function [OBJ,V,DELTA,OUTP] = mykalman(THIS,INP,OUTP,OPT,varargin)
% mykalman  [Not a public function] Kalman filter.
%
% Backed IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

% This Kalman filter handles the following special cases:
% * non-stationary initial conditions (treated as fixed numbers);
% * measurement parameters concentrated out of likelihood;
% * time-varying std deviations;
% * exclusion of some of the measurement variables from likelihood;
% * exclusion of some of the periods from likelihood;
% * k-step-ahead predictions;
% * rollbacks;
% * tunes on the mean of shocks, combined anticipated and unanticipated;
% * missing observations entered as NaNs;
% * infinite std devs of measurement shocks equivalent to missing obs.

NARGOUT = nargout;
ny = length(THIS.solutionid{1});
nx = size(THIS.solution{1},1);
nb = size(THIS.solution{1},2);
nf = nx - nb;
ne = length(THIS.solutionid{3});
neqtn = length(THIS.eqtn);
nalt = size(THIS.Assign,3);
ndata = size(INP,3);

% TODO: Output data for the updating step.
%**************************************************************************

s = struct();

s.isnonlin = OPT.nonlinearise > 0 && any(THIS.nonlin);
s.ahead = OPT.ahead;

realsmall = getrealsmall();

% Out-of-lik params cannot be used with ~opt.dtrends.
npout = length(OPT.outoflik);

% Source data are always numeric arrays without initial condition.
% Add pre-sample initial condition.
INP = [nan(ny,1,ndata),INP];
nper = size(INP,2);

% Struct with currently processed information. Initialise the invariant
% fields.
s.ny = ny;
s.nx = nx;
s.nb = nb;
s.nf = nf;
s.ne = ne;
s.nalt = nalt;
s.nper = nper;
s.npout = npout;

% Add presample to objective function range and time trend.
s.objrange = [false,OPT.objrange];
s.ttrend = [nan,OPT.ttrend];

% Do not adjust the option `'lastsmooth='` -- see comments in `loglikopt`.
s.lastsmooth = OPT.lastsmooth;

% Tunes on shock means; model solution is expanded within `mypreloglik`.
tune = OPT.tune;
s.istune = ~isempty(tune) && any(tune(:) ~= 0);
if s.istune
    % Add pre-sample.
    ntune = size(tune,3);
    tune = [zeros(ne,1,ntune),tune];
end

% Total number of cycles.
nloop = max(ndata,nalt);

% Pre-allocate first two outpu args.
OBJ = nan(1,nloop,OPT.precision);
V = nan(1,nloop,OPT.precision);
DELTA = nan(1,npout,nloop,OPT.precision);

% Pre-allocate output data.
dorequestoutp();
if NARGOUT >= 4
    dopreallocoutp();
end

% Logical indices of shocks occuring in measurement and transition
% equations.
s = xxshocktypes(s,THIS);

% Prepare struct and options for non-linear simulations (prediction step).
if s.isnonlin
    dopreparenonlin();
end

% MAIN LOOP
% ==========

if OPT.progress
    progress = progressbar('IRIS model.kalman progress');
end

for iloop = 1 : nloop
    
    % Next model solution.
    if iloop <= nalt
        [T,R,k,s.Z,s.H,s.d,s.U] = mysspace(THIS,iloop,true);
        s.abseigval = abs(ordeig(T(nf+1:end,:)));
        s.nunit = sum(abs(s.abseigval - 1) <= realsmall);
        s.Tf = T(1:nf,:);
        s.Ta = T(nf+1:end,:);
        % Keep forward expansion for computing the effect of tunes on shock means.
        % Cut off the expansion within subfunctions.
        s.Rf = R(1:nf,:);
        s.Ra = R(nf+1:end,:);
        s.Zt = s.Z.';
        if OPT.deviation
            s.ka = [];
            s.kf = [];
            s.d = [];
        else
            s.d = THIS.solution{6}(:,:,iloop);
            s.k = THIS.solution{3}(:,1,iloop);
            s.kf = k(1:nf,:);
            s.ka = k(nf+1:end,:);
        end
        
        % TIME-VARYING STDCORR
        % =====================
        % Combine currently assigned `stdcorr` in the model object with the
        % user-supplied time-vaying `stdcorr`.
        thisstdcorr = THIS.stdcorr(1,:,iloop).';
        s.stdcorr = metaobj.mycombinestdcorr(thisstdcorr,OPT.stdcorr,nper-1);

        % Add presample, which will be used to initialise the Kalman
        % filter.
        s.stdcorr = [thisstdcorr,s.stdcorr];
        
        % Create covariance matrix from stdcorr vector.
        s.Omg = covfun.stdcorr2cov(s.stdcorr,ne);
        
        % Create reduced form covariance matrices `Sa` and `Sy`, and find
        % measurement variables with infinite measurement shocks, `syinf`.
        s = xxomg2sasy(s);
        
        % Free memory.
        s.stdcorr = [];
        
        % DETERMINISTIC TRENDS
        % =====================
        % y(t) - D(t) - X(t)*delta = Z*a(t) + H*e(t).
        if npout > 0 || OPT.dtrends
            [s.D,s.X] = mydtrends(THIS,s.ttrend,OPT.outoflik,iloop);
        else
            s.D = [];
            s.X = zeros(ny,0,nper);
        end
        
        % INITIAL DISTRIBUTION
        % =====================
        doinitcond();
        
    end
    
    % Next tunes on the means of the shocks.
    if s.istune && iloop <= ntune
        s.tune = tune(:,:,iloop);
    end
    % Add the effect of the tunes to the constant vector; recompute the
    % effect whenever the tunes have changed or the model solution has
    % changed or both.
    %
    % The std dev of the tuned shocks remain unchanged and
    % hence the filtered shocks can differ from its tunes (unless the user
    % specifies zero std dev).
    if s.istune && (iloop <= ntune || iloop <= nalt)
        [s.d,s.ka,s.kf] = xxshocktunes(s,OPT);
    end
    
    % Next data.
    if iloop <= ndata
        s.data = INP(:,:,iloop);
    end
    % Subtract fixed deterministic trends from measurement variables.
    s.y1 = s.data;
    if ~isempty(s.D)
        s.y1 = s.y1 - s.D;
    end
    
    % Make measurement variables with `Inf` measurement shocks look like
    % missing. The `Inf` measurement shocks are detected in `sub_omg2sasy`.
    s.y1(s.syinf) = NaN;
    
    % Index of available observations.
    s.yindex = ~isnan(s.y1);
    s.jyeq = [false,all(s.yindex(:,2:end) == s.yindex(:,1:end-1),1)];
    
    % NUMBER OF INITIAL CONDITIONS TO BE OPTIMISED
    % =============================================
    if isequal(OPT.initcond,'optimal')
        % All initial conditions, including stationary variables, will be
        % optimised.
        s.ninit = nb;
    elseif iscell(OPT.initcond)
        % No initial condition will be optimised if it is supplied by the user.
        s.ninit = 0;
    else
        % Estimate diffuse initial conditions only if there is at least one
        % diffuse measurement variable with at least one observation.
        z = s.Z(any(s.yindex,2),1:s.nunit);
        if any(any(abs(z) > realsmall))
            s.ninit = s.nunit;
        else
            s.ninit = 0;
        end
    end
    
    % Run prediction error decomposition and evaluate user-requested
    % objective function.
    [OBJ(iloop),s] = xxped(THIS,s,OPT,iloop);
    
    % Correct prediction step for estimated initial conditions and dtrends
    % parameters.
    if s.storeped
        s = xxcorrectpe(s);
    end
    % Free memory.
    s.M = [];
    
    if s.storepredict && (s.ninit > 0 || npout > 0)
        s = xxcorrect(s);
    end
    % Free memory.
    s.Q = [];
    
    % Calculate prediction step for fwl variables.
    if s.retpred || s.retsmooth
        s.isxfmse = s.retpredstd || s.retpredmse ...
            || s.retsmoothstd || s.retsmoothmse;
        s = xxfilterxfmean(s);
    end
    
    % K-step-ahead predictions.
    if s.ahead > 1 && s.storepredict
        s = xxahead(s);
    end
    
    % Non-essential reporting MSEs and stds in prediction step.
    if s.retpredstd || s.retpredmse ...
            || s.retsmoothstd || s.retsmoothmse
        s = xxpredmse(s);
    end
    
    % Smoothing step for all variables.
    if s.retsmooth
        s.lastobs = find(any(s.yindex,1),1,'last');
        if s.retsmoothstd || s.retsmoothmse
            s = xxsmoothmse(s);
        end    
        s = xxsmoothmean(s);
    end
    
    % Columns in `pe` to be filled.
    if s.ahead > 1
        predcols = 1 : s.ahead;
    else
        predcols = iloop;
    end
    
    if s.retpred
        % Add prediction data to output data handle.
        doretpred();
    end
    
    if s.retsmooth
        % Add smoothing data to output data handle.
        doretsmooth();
    end
    
    if s.storeped
        V(iloop) = s.V;
        DELTA(1,:,iloop) = s.delta;
        if s.retpe
            OUTP.pe(:,:,predcols) = s.pe;
        end
        if s.retf
            OUTP.F(:,:,:,iloop) = s.F*s.V;
        end
        if s.retpdelta
            OUTP.Pdelta(:,:,iloop) = s.Pdelta*s.V;
        end
    end
    
    if OPT.progress
        update(progress,iloop/nloop);
    end
    
end

% Nested functions.

%**************************************************************************
    function dorequestoutp()
        % dooutprequest  Check for requested output data.
        s.retpe = NARGOUT >= 3 && isfield(OUTP,'pe');
        s.retf = NARGOUT >= 3 && isfield(OUTP,'F');
        s.retpdelta = NARGOUT >= 3 && isfield(OUTP,'Pdelta');
        s.retpredmean = NARGOUT >= 3 && isfield(OUTP,'predmean');
        s.retpredmse = NARGOUT >= 3 && isfield(OUTP,'predmse');
        s.retpredstd = NARGOUT >= 3 && isfield(OUTP,'predstd');
        s.retsmoothmean = NARGOUT >= 3 && isfield(OUTP,'smoothmean');
        s.retsmoothmse = NARGOUT >= 3 && isfield(OUTP,'smoothmse');
        s.retsmoothstd = NARGOUT >= 3 && isfield(OUTP,'smoothstd');
        s.retpred = s.retpredmean || s.retpredstd || s.retpredmse;
        s.retsmooth = s.retsmoothmean || s.retsmoothstd || s.retsmoothmse;
        s.storeped = NARGOUT >= 3 && ( ...
            s.retpdelta || s.retpe ...
            || s.retpred || s.retsmooth);
        s.storepredict = s.ahead > 1 || s.retpred || s.retsmooth;
    end

%**************************************************************************
    function dopreallocoutp()
        % dopreallocoutp  Preallocate output data.
        % This function is called only when `nargout` >= 3.
        if s.retpdelta
            OUTP.Pdelta = nan(npout,npout,nloop,OPT.precision);
        end
        npred = max(nloop,s.ahead);
        if s.retpe
            % `DATA.pe` and `DATA.predmean` may include k-step-ahead predictions.
            OUTP.pe = nan(ny,nper,npred,OPT.precision);
        end
        if s.retf
            OUTP.F = nan(ny,ny,nper,nloop,OPT.precision);
        end
        if s.retpredmean
            % `DATA.pe` and `DATA.predmean` may include k-step-ahead predictions.
            OUTP.predmean = myhdatainit(THIS,nper,npred,OPT.precision);
        end
        if s.retpredstd
            OUTP.predstd = myhdatainit(THIS,nper,nloop,OPT.precision);
        end
        if s.retpredmse
            OUTP.predmse = nan(nb,nb,nper,nloop,OPT.precision);
        end
        if s.retsmoothmean
            OUTP.smoothmean = myhdatainit(THIS,nper,nloop,OPT.precision);
        end
        if s.retsmoothstd
            OUTP.smoothstd = myhdatainit(THIS,nper,nloop,OPT.precision);
        end
        if s.retsmoothmse
            OUTP.smoothmse = nan(nb,nb,nper,nloop,OPT.precision);
        end        
    end
% dopreallocoutp().

%**************************************************************************
    function doretpred()
        % Return pred mean.
        % Note that s.y0, s.f0 and s.a0 include k-sted-ahead predictions if
        % ahead > 1.
        if s.retpredmean
            yy = s.y0;
            % Convert `alpha` predictions to `xb` predictions. The
            % `a0` may contain k-step-ahead predictions in 3rd dimension.
            bb = s.a0;
            for ii = 1 : size(bb,3)
                bb(:,:,ii) = s.U*bb(:,:,ii);
            end
            xx = [s.f0;bb];
            % Shock predictions are always zeros.
            ee = zeros(ne,nper,s.ahead);
            % Set predictions for the pre-sample period to `NaN`.
            yy(:,1,:) = NaN;
            xx(:,1,:) = NaN;
            ee(:,1,:) = NaN;
            % Add fixed deterministic trends back to measurement vars.
            if ~isempty(s.D)
                % We need to use `bsxfun` to get s.D expanded along 3rd dimension if
                % `ahead=` is greater than 1. This is equivalent to the older version:
                %     yy = yy + s.D(:,:,ones(1,ahead));
                yy = bsxfun(@plus,yy,s.D);
            end
            % Add shock tunes to shocks.
            if s.istune
                % This line is equivalent to the older version:
                %     ee = ee + s.tune(:,:,ones(1,ahead));
                ee = bsxfun(@plus,ee,s.tune);
            end
            % Do not use lags in the prediction output data.
            myhdataassign(THIS,OUTP.predmean,predcols, ...
                yy, ...
                xx, ...
                ee, ...
                '-nopresample');
        end
        
        % Return pred std.
        if s.retpredstd
            % Do not use lags in the prediction output data.
            myhdataassign(THIS,OUTP.predstd,iloop, ...
                s.Dy0*s.V, ...
                [s.Df0;s.Db0]*s.V, ...
                s.De0*s.V, ...
                '-nopresample','-std');
        end

        if s.retpredmse
            OUTP.predmse(:,:,:,iloop) = s.Pb0*s.V;
        end

    end
% doretpred().

%**************************************************************************
    function doretsmooth()
        % Return smooth mean.
        if s.retsmoothmean
            yy = s.y2;
            xx = [s.f2;s.U*s.a2(:,:,1)];
            ee = s.e2;
            yy(:,1:s.lastsmooth) = NaN;
            xx(:,1:s.lastsmooth-1) = NaN;
            xx(1:nf,s.lastsmooth) = NaN;
            ee(:,1:s.lastsmooth) = NaN;
            % Add deterministic trends to measurement vars.
            if ~isempty(s.D)
                yy = yy + s.D;
            end
            % Add shock tunes to shocks.
            if s.istune
                ee = ee + s.tune;
            end
            myhdataassign(THIS,OUTP.smoothmean,iloop,yy,xx,ee);
        end
        
        % Return smooth std.
        if s.retsmoothstd
            s.Dy2(:,1:s.lastsmooth) = NaN;
            s.Df2(:,1:s.lastsmooth) = NaN;
            s.Db2(:,1:s.lastsmooth-1) = NaN;
            myhdataassign(THIS,OUTP.smoothstd,iloop, ...
                s.Dy2*s.V, ...
                [s.Df2;s.Db2]*s.V, ...
                [], ...
                '-std');
        end
        
        % Return smooth MSE for `xb`.
        if s.retsmoothmse
            s.Pb2(:,:,1:s.lastsmooth-1) = NaN;
            OUTP.smoothmse(:,:,:,iloop) = s.Pb2*s.V;
        end
    end
% doretsmooth().

%**************************************************************************
    function doinitcond()
        % doinitcond  Set up initial condition.
        % Set up initial condition for the mean and MSE matrix.
        s.Painit = zeros(nb);
        tmpstable = [false(1,s.nunit),true(1,nb-s.nunit)];
        % Initialise MSE matrix.
        if iscell(OPT.initcond)
            % User-supplied initial condition.
            s.Painit = OPT.initcond{2}(:,:,1,min(end,iloop));
        elseif nb > s.nunit && isequal(OPT.initcond,'stochastic')
            % R matrix with rows corresponding to stable alpha and columns
            % corresponding to transition shocks.
            RR = s.Ra(:,1:ne);
            RR = RR(tmpstable,s.tshocks);
            % Reduced form covariance corresponding to stable alpha. Use the structural
            % shock covariance sub-matrix corresponding to transition shocks only in
            % the pre-sample period.
            Sa = RR*s.Omg(s.tshocks,s.tshocks,1)*RR.';
            % Compute asymptotic initial condition.
            if sum(tmpstable) == 1
                Pa0stable = Sa / (1 - s.Ta(tmpstable,tmpstable)^2);
            else
                Pa0stable = ...
                    covfun.lyapunov(s.Ta(tmpstable,tmpstable),Sa);
                Pa0stable = (Pa0stable + Pa0stable.')/2;
            end
            s.Painit(tmpstable,tmpstable) = Pa0stable;
        end
        % Initialise mean.
        s.ainit = zeros(nb,1);
        if iscell(OPT.initcond)
            % User-supplied initial condition.
            s.ainit = OPT.initcond{1}(:,1,min(end,iloop));
        elseif ~isempty(s.ka)
            % Asymptotic initial condition.
            s.ainit(tmpstable,1) = ...
                (eye(nb-s.nunit) - s.Ta(tmpstable,tmpstable)) ...
                \ s.ka(tmpstable,1);
        end
    end
% doinitcond().

%**************************************************************************
    function dopreparenonlin()
        s2.isplan = false;
        s2.istune = false;
        s2.isnonlinplan = false;
        s2.isnonlinopt = true;
        s2.isnonlin = true;
        s2.lastea = 0;
        s2.lasteu = 0;
        s2.qanchors = false(neqtn,OPT.nonlinearise);
        s2.qanchors(THIS.nonlin,:) = true;
        s2.label = myget(THIS,'canbenonlinearised');
        s2.lastendog = 0;
        s2.lastendogu = 0;
        s2.lastexog = 0;
        s2.lastnonlin = OPT.nonlinearise;
        s2.tplusk = s2.lastnonlin - 1;
        s2.progress = [];
        s2.a0 = [];
        s2.e = zeros(ne,1);%OPT.nonlinearise);
        s2.ytune = [];
        s2.xtune = [];
        s2.W = [];
        s2.zerothsegment = 0;
        s.s2 = s2;
        s.simulateopt = passvalopt('model.simulate',OPT.simulate{:});
    end
% dopreparenonlin().

end

% Subfunctions.

%**************************************************************************
function s = xxcorrectpe(s)
% xxcorrectpe  Correct prediction errors for the estimated init and delta.
est = [s.delta(:);s.init(:)];
nper = size(s.y1,2);
for t = 2 : nper
    j = s.yindex(:,t);
    s.pe(j,t) = s.pe(j,t) - s.M(j,:,t)*est;
end
end
% xxcorrectpe().

%**************************************************************************
function s = xxcorrect(s)
% xxcorrect  Correct alpha, y, and F\pe for estimated initial conditions, `init`, and
% out-of-lik parameters, and deltas.
[ny,nper] = size(s.y1);
init = s.init(:);
delta = s.delta(:);
est = [delta;init];
% Store the effect of out-of-lik parameters, `delta`, on measurement
% variables because we need to correct k-step ahead predictions and
% smoothed estimates. The effect of diffuse init conditions, `init`,
% will have been already accounted for in the estimates of `alpha`.
s.ydelta = nan(ny,nper);
s.a0(:,1) = s.a0(:,1) + s.Q(:,:,1)*est;
for t = 2 : nper
    correct = s.Q(:,:,t)*est; % s.Q1(:,:,t)*delta + s.Q2(:,:,t)*init;
    s.a0(:,t) = s.a0(:,t) + correct;
    s.ydelta(:,t) = s.X(:,:,t)*delta;
end
s.y0 = s.Z*s.a0 + s.ydelta;
if ~isempty(s.d)
    % We need to use `bsxfun` because `s.d` can be either a column vector (if
    % there are no tunes on the shock means) or a full array (if there are
    % tunes on shock means). This is equivalent to the older version:
    %     if ~s.istune
    %         s.y0 = s.y0 + s.d(:,ones(1,nper));
    %     else
    %         s.y0 = s.y0 + s.d;
    %     end
    s.y0 = bsxfun(@plus,s.y0,s.d);
end
end
% xxcorrect().

%**************************************************************************
function s = xxahead(s)
% xxahead  K-step ahead predictions and prediction errors for K>2 when
% requested by caller. This function must be called after correction for
% diffuse initial conditions and/or out-of-lik params has been made.
s.a0 = cat(3,s.a0,nan([size(s.a0),s.ahead-1]));
s.pe = cat(3,s.pe,nan([size(s.pe),s.ahead-1]));
s.y0 = cat(3,s.y0,nan([size(s.y0),s.ahead-1]));
if s.retpred
    % `s.f0` exists and needs to be calculated only if `pred` data are
    % requested.
    s.f0 = cat(3,s.f0,nan([size(s.f0),s.ahead-1]));
end
nper = size(s.a0,2);
for k = 2 : min(s.ahead,nper-1);
    t = 1+k : nper;
    tmpindex = ones(1,numel(t));
    s.a0(:,t,k) = s.Ta*s.a0(:,t-1,k-1);
    if ~isempty(s.ka)
        if ~s.istune
            s.a0(:,t,k) = s.a0(:,t,k) + s.ka(:,tmpindex);
        else
            s.a0(:,t,k) = s.a0(:,t,k) + s.ka(:,t);
        end
    end
    s.y0(:,t,k) = s.Z*s.a0(:,t,k);
    if ~isempty(s.d)
        if ~s.istune
            s.y0(:,t,k) = s.y0(:,t,k) + s.d(:,tmpindex);
        else
            s.y0(:,t,k) = s.y0(:,t,k) + s.d(:,t);
        end
    end
    if s.retpred
        s.f0(:,t,k) = s.Tf*s.a0(:,t-1,k-1);
        if ~isempty(s.kf)
            if ~s.istune
                s.f0(:,t,k) = s.f0(:,t,k) + s.kf(:,tmpindex);
            else
                s.f0(:,t,k) = s.f0(:,t,k) + s.kf(:,t);
            end
        end
    end
end
if s.npout > 0
    s.y0(:,:,2:end) = s.y0(:,:,2:end) + s.ydelta(:,:,ones(1,s.ahead-1));
end
s.pe(:,:,2:end) = s.y1(:,:,ones(1,s.ahead-1)) - s.y0(:,:,2:end);
end
% xahead().

%**************************************************************************
function [OBJ,S] = xxped(THIS,S,OPT,ILOOP)
% xxped  Prediction step and evaluation of objective function.

% Size of system matrices.
ny = S.ny;
nb = S.nb;
nf = S.nf;
ne = S.ne;
nper = size(S.y1,2);
npout = S.npout;
ninit = S.ninit;
lastomg = size(S.Omg,3);

Ta = S.Ta;
Tat = S.Ta.';
ka = S.ka;
d = S.d;
jy = false(ny,1);
Z = S.Z(jy,:);
X = S.X;

y1 = S.y1;
yindex = S.yindex;
a = S.ainit;
P = S.Painit;
G = zeros(nb,0);
K = zeros(nb,0);
pe = zeros(0,1);
ZP = zeros(0,nb);

% Initialise objective function components.
sumpeFipe = 0;
sumlogdetF = 0;
% Effect of outofliks and fixed init states on a(t).
Q1 = zeros(nb,npout);
Q2 = eye(nb,ninit);
% Effect of outofliks and fixed init states on pe(t).
M1 = zeros(0,npout);
M2 = zeros(0,ninit);

% Initialise flags.
anypout = npout > 0;
anyinit = ninit > 0;
anyestimate = anypout || anyinit;
objrange = S.objrange;

% Initialise sum terms used in out-of-lik estimation.
if anyestimate
    sumMtFiM = 0;
    sumMtFipe = 0;
else
    sumMtFiM = [];
    sumMtFipe = [];
end

% Initialise matrices that are to be stored.
if S.storeped
    S.F = nan(ny,ny,nper);
    S.Fd = nan(1,nper);
    S.pe = nan(ny,nper);
    S.M = nan(ny,npout+ninit,nper);
    if S.storepredict
        S.a0 = nan(nb,nper);
        S.a0(:,1) = S.ainit;
        S.Pa0 = nan(nb,nb,nper);
        S.Pa0(:,:,1) = S.Painit;
        S.Pa1 = S.Pa0;
        S.De0 = nan(ne,nper);
        S.y0 = nan(ny,nper);
        S.K = nan(nb,ny,nper);
        S.Q = zeros(nb,npout+ninit,nper);
        S.Q(:,npout+1:end,1) = Q2;
        if S.retsmooth
            S.L = nan(nb,nb,nper);
            S.L(:,:,1) = Ta;
        end
    end
end

% Free memory.
S.Painit = [];
S.ainit = [];

% Main loop.
for t = 2 : nper
    
    % Effect of outofliks on `a(t)`. This step must be made before
    % updating `jy` because we use `G` from the previous period.
    if anypout
        Q1 = Ta*Q1 - G*M1(jy,:);
    end
    
    % Effect of fixed init states on `a(t)`. This step must be made
    % before updating `jy` because we use `G` from the previous
    % period.
    if anyinit
        Q2 = Ta*Q2 - G*M2(jy,:);
    end
    
    if ~S.isnonlin
        % Prediction `a(t|t-1)` based on `a(t-1|t-2)`, prediction error
        % `pe(t-1)`, and the Kalman gain `G(t-1)`.
        a = Ta*a + G*pe;
        % Adjust the prediction step for the constant vector.
        if ~isempty(ka)
            if ~S.istune
                a = a + ka;
            else
                a = a + ka(:,t);
            end
        end
    else
        % Run non-linear simulation to produce the mean prediction.
        dononlinpredict();
    end
    
    % Update reduced-form shock covariance if the structural shock
    % covariance has changed from t-1. We need to update `Omg` before
    % computing the index of available observations, `jy`, because
    % possible infinite variances in `Omg` are translated into missing
    % observations.
    tomg = min(t,lastomg);
    Omg = S.Omg(:,:,tomg);
    Sa = S.Sa(:,:,tomg);
    Sy = S.Sy(:,:,tomg);
    
    % Index of available observations.
    jy = yindex(:,t);
    
    % Amongst availabe observations, these are not excluded.
    include = ~OPT.exclude(jy);
    anyexclude = any(~include);
    
    % Prediction MSE P(t|t-1) based on P(t-1|t-2) and the predictive Kalman
    % gain G(t-1).
    P = (Ta*P - G*ZP)*Tat + Sa;
    % Make sure P is symmetric and does not explode over time.
    P = (P+P.')/2;
    % Time-varying Z matrix.
    Z = S.Z(jy,:);
    % Prediction MSE for observables F(t).
    ZP = Z*P;
    PZt = ZP.';
    F = Z*PZt + Sy(jy,jy);
    K = PZt/F;
    G = Ta*K;
    
    if anyexclude
        % `Fx` excludes variables requested by the user from computing the
        % likelihood function.
        Fx = F(include,include);
        Fd = log(det(Fx));
    else
        Fd = log(det(F));
    end
    
    % Prediction error `pe(t)`. The size of pe changes with `jy`.
    pe = y1(jy,t) - Z*a;
    % Adjust the predicition error for the constant vector.
    if ~isempty(d)
        if ~S.istune
            pe = pe - d(jy,:);
        else
            pe = pe - d(jy,t);
        end
    end
    
    if anyestimate
        % Update effect of outofliks and fixed init states on pe(t).
        M1 = S.Z*Q1 + X(:,:,t);
        M2 = S.Z*Q2;
        M = [M1,M2];
        Mt = M.';
    end
    
    if objrange(t)
        iy = ~OPT.exclude & jy;
        % Compute components of objective functions if user wants to
        % include this period of observations.
        if OPT.objective == 1
            % Likelihood function.
            if anyexclude
                % If some observations excluded, we must divide pred
                % errors by a submatrix of F.
                pex = pe(include);
                pext = pex.';
                sumpeFipe = sumpeFipe + pext*(Fx\pex);
            else
                pet = pe.';
                sumpeFipe = sumpeFipe + (pet/F)*pe;
            end
            sumlogdetF = sumlogdetF + Fd;
            % Estimation of pout and delta.
            if anyestimate
                if anyexclude
                    MtFi = Mt(:,iy)/Fx;
                    sumMtFipe = sumMtFipe + MtFi*pex;
                else
                    MtFi = Mt(:,iy)/F;
                    sumMtFipe = sumMtFipe + MtFi*pe;
                end
                sumMtFiM = sumMtFiM + MtFi*M(iy,:);
            end
        elseif OPT.objective == 2
            % Weighted sum of prediction errors.
            pex = pe(include);
            pext = pe.';
            W = OPT.weighting(iy,iy);
            sumpeFipe = sumpeFipe + pext*W*pex;
            % Estimation of pout and delta.
            if anyestimate
                MtFi = Mt(:,iy)*W;
                sumMtFipe = sumMtFipe + MtFi*pex;
                sumMtFiM = sumMtFiM + MtFi*M(iy,:);
            end
        end
    end
    
    if S.storeped
        if S.storepredict
            % Compute MSE for all measurement variables, not only for
            % the currently observed ones when predict data are
            % returned.
            S.F(:,:,t) = S.Z*P*S.Z.' + Sy;
        else
            S.F(jy,jy,t) = F;
        end
        S.pe(jy,t) = pe;
        if anyestimate
            S.M(:,:,t) = M;
        end
        if S.storepredict
            dostorepredict();
        end
    end
    
end
% End of time loop.

% Total number of observations included in loglik.
nobs = sum(sum(...
    yindex(:,2:end) & ~OPT.exclude(:,ones(1,nper-1))...
    ));

% Evaluate likelihood, variance factor, out-of-lik parameters.
[OBJ,V,est,Pest] = ...
    model.myoutoflik(sumlogdetF,sumpeFipe,sumMtFiM,sumMtFipe,nobs,OPT);
S.delta = est(1:npout);
S.Pdelta = Pest(1:npout,1:npout);
S.init = est(npout+1:end);
S.V = V;

% Free memory.
S.Sa = [];
S.Sy = [];

% Nested functions.

%**************************************************************************
    function dostorepredict()
        % Store all necessary information on this prediction step.
        S.a0(:,t) = a;
        S.Pa0(:,:,t) = P;
        S.Pa1(:,:,t) = P - K*ZP;
        S.De0(:,t) = diag(Omg);
        S.y0(:,t) = S.Z*a;
        if ~isempty(d)
            if ~S.istune
                S.y0(:,t) = S.y0(:,t) + d;
            else
                S.y0(:,t) = S.y0(:,t) + d(:,t);
            end
        end
        S.K(:,jy,t) = K;
        S.Q(:,:,t) = [Q1,Q2];
        if S.retsmooth
            S.L(:,:,t) = Ta - G*Z;
        end
    end
% dostorepredict().

%**************************************************************************
    function dononlinpredict()
        a1 = a + K*pe;
        S.s2.a0 = a1;
        S.s2.zerothsegment = t-2;
        s2 = mysimulate(THIS,S.s2,S.simulateopt,ILOOP);
        a = s2.w(nf+1:end,1);
    end

end
% xxped().

%**************************************************************************
function s = xxfilterxfmean(s)
% xxfilterxfmean  Point prediction step for fwl transition variables. The
% MSE matrices are computed when needed in `xxsmoothmse` on the spot.
nf = s.nf;
nper = s.nper;
% Pre-allocate state vectors.
s.f0 = nan(nf,nper);
if nf == 0
    return
end
for t = 2 : nper
    j1 = s.yindex(:,t-1);
    s.f0(:,t,1) = s.Tf*(s.a0(:,t-1) + s.K(:,j1,t-1)*s.pe(j1,t-1));
    if ~isempty(s.kf)
        if ~s.istune
            s.f0(:,t,1) = s.f0(:,t,1) + s.kf;
        else
            s.f0(:,t,1) = s.f0(:,t,1) + s.kf(:,t);
        end
    end
end
end
% xxfilterxfmean().

%**************************************************************************
function s = xxpredmse(s)
s.Pb0 = xxpa2pb(s.U,s.Pa0);
if s.retpredstd || s.retsmoothstd
    s.Dy0 = nan(s.ny,s.nper);
    s.Df0 = nan(s.nf,s.nper);
    s.Db0 = nan(s.nb,s.nper);
    for t = 2 : s.nper
        s.Dy0(:,t) = diag(s.F(:,:,t));        
        s.Db0(:,t) = diag(s.Pb0(:,:,t));
    end
    % Compute `Df0`.
    s = xxfilterxfmse(s,2:s.nper);
end
end
% xxpredmse().

%**************************************************************************
function s = xxsmoothmse(s)
% xxsmoothmse  Smoother for MSE matrices of all variables.
ny = s.ny;
nb = s.nb;
nf = s.nf;
nper = s.nper;
lastsmooth = s.lastsmooth;
U = s.U;
% Pre-allocation.
if s.retsmoothmse
    s.Pb2 = nan(nb,nb,nper);
end
s.Db2 = nan(nb,nper); % Diagonal of Pb2.
s.Df2 = nan(nf,nper); % Diagonal of Pf2.
s.Dy2 = nan(ny,nper); % Diagonal of Py2.
% Index of available observations.
lastobs = find(any(s.yindex,1),1,'last');
if lastobs < nper
    s.Pb2(:,:,lastobs+1:nper) = s.Pb0(:,:,lastobs+1:nper);
    s.Dy2(:,lastobs+1:nper) = s.Dy0(:,lastobs+1:nper);
    s.Df2(:,lastobs+1:nper) = s.Df0(:,lastobs+1:nper);
    s.Db2(:,lastobs+1:nper) = s.Db0(:,lastobs+1:nper);
end
N = 0;
for t = lastobs : -1 : lastsmooth
    j = s.yindex(:,t);
    N = (s.Z(j,:).'/s.F(j,j,t))*s.Z(j,:) + s.L(:,:,t).'*N*s.L(:,:,t);
    Pa0NPa0 = s.Pa0(:,:,t)*N*s.Pa0(:,:,t);
    Pa2 = s.Pa0(:,:,t) - Pa0NPa0;
    Pa2 = (Pa2+Pa2.')/2;
    Pb2 = xxpa2pb(U,Pa2);
    if s.retsmoothmse
        s.Pb2(:,:,t) = Pb2;
    end
    s.Db2(:,t) = diag(Pb2);
    if nf > 0 && t > lastsmooth
        % Do not run smoother on `xf` in the very first period.
        [~,Pf0,Pfa0] = xxfilterxfmse(s,t);
        % Fwl transition variables.
        Pfa0N = Pfa0*N;
        Pf2 = Pf0 - Pfa0N*Pfa0.';
        % Pfa2 = s.Pfa0(:,:,t) - Pfa0N*s.Pa0(:,:,t);
        Pf2 = (Pf2+Pf2.')/2;
        s.Df2(:,t) = diag(Pf2);
    end
    if ny > 0
        % Measurement variables.
        Py2 = s.F(:,:,t) - s.Z*Pa0NPa0*s.Z.';
        Py2 = (Py2+Py2.')/2;
        Py2(j,:) = 0;
        Py2(:,j) = 0;
        s.Dy2(:,t) = diag(Py2);
    end   
end
end
% xxsmoothmse().

%**************************************************************************
function [s,Pf0,Pfa0] = xxfilterxfmse(s,time)
% xxfilterxfmse  MSE prediction step for fwl variables.
ne = s.ne;
nf = s.nf;
if nf == 0
    return
end
lastomg = size(s.Omg,3);
tshocks = s.tshocks;
% Cut off forward expansion.
Rf = s.Rf(:,1:ne);
Ra = s.Ra(:,1:ne);
for t = time
    tomg = min(t,lastomg);
    TfPa1 = s.Tf*s.Pa1(:,:,t-1);
    RfOmg = Rf(:,tshocks)*s.Omg(tshocks,:,tomg);
    Pf0 = TfPa1*s.Tf.' + RfOmg*Rf.';
    Pf0 = (Pf0 + Pf0.')/2;
    Pfa0 = TfPa1*s.Ta.' + RfOmg*Ra.';
    if s.retpredstd || s.retsmoothstd
        s.Df0(:,t) = diag(Pf0);
    end
end
end
% xxfilterxfmse().

%**************************************************************************
function s = xxsmoothmean(s)
% xxsmoothmean  Kalman smoother for point estimates of all variables.
nx = s.nx;
nb = s.nb;
nf = nx - nb;
ne = s.ne;
nper = size(s.y1,2);
npout = s.npout;
tshocks = s.tshocks;
mshocks = s.mshocks;
lastomg = size(s.Omg,3);
% Cut off forward expansion.
Rf = s.Rf(:,1:ne);
Ra = s.Ra(:,1:ne);
% Pre-allocation. Re-use first page of prediction data. Prediction data
% can have multiple pages if ahead > 1.
s.a2 = s.a0(:,:,1);
s.f2 = s.f0(:,:,1);
s.e2 = zeros(ne,nper);
s.y2 = s.y1(:,:,1);
% Last available observation.
% Do not run the smoother beyond that point.
lastobs = max(0,find(any(s.yindex,1),1,'last'));
s.y2(:,lastobs+1:end) = s.y0(:,lastobs+1:end,1);
r = zeros(nb,1);
for t = lastobs : -1 : s.lastsmooth
    j = s.yindex(:,t);
    Fipe = s.F(j,j,t) \ s.pe(j,t,1);
    tomg = min(t,lastomg);
    if any(j)
        G = s.Ta*s.K(:,j,t);
        HOmg = s.H(j,mshocks)*s.Omg(mshocks,:,tomg);
        s.e2(:,t) = s.e2(:,t) ...
            + HOmg.'*(Fipe - G.'*r);
    end
    r = s.Zt(:,j)*Fipe + s.L(:,:,t).'*r;
    s.a2(:,t) = s.a0(:,t,1) + s.Pa0(:,:,t)*r;
    RaOmg = Ra(:,tshocks)*s.Omg(tshocks,:,tomg);
    s.e2(:,t) = s.e2(:,t) + RaOmg.'*r;
    % Back out NaN observables.
    if any(~j)
        s.y2(~j,t) = s.Z(~j,:)*s.a2(:,t) + s.H(~j,:)*s.e2(:,t);
        if ~isempty(s.d)
            % Correct the estimates of NaN observations for the
            % deterministic trends.
            if ~s.istune
                s.y2(~j,t) = s.y2(~j,t) + s.d(~j,:);
            else
                s.y2(~j,t) = s.y2(~j,t) + s.d(~j,t);
            end
        end
        if npout > 0
            % Correct the estimates of NaN observations for the effect
            % of estimated out-of-lik parameters.
            s.y2(~j,t) = s.y2(~j,t) + s.ydelta(~j,t);
        end
    end
end
% Fwl transition variables.
if nf > 0
    s.f2(:,2:nper) = s.Tf*s.a2(:,1:nper-1) + Rf*s.e2(:,2:nper);
    if ~isempty(s.kf)
        if ~s.istune
            s.f2(:,2:nper) = s.f2(:,2:nper) + s.kf(:,ones(1,nper-1));
        else
            s.f2(:,2:nper) = s.f2(:,2:nper) + s.kf(:,2:nper);
        end
    end
end
end
% xxsmoothmean().

%**************************************************************************
function [d,ka,kf] = xxshocktunes(s,opt)
% xxshocktunes  Add tunes on shock means to constant terms.
ny = s.ny;
nb = s.nb;
nf = s.nf;
ne = s.ne;
nper = s.nper;
if ne == 0
    return
end
H = s.H;
Rf = s.Rf;
Ra = s.Ra;
if opt.deviation
    d = zeros(ny,nper);
    ka = zeros(nb,nper);
    kf = zeros(nf,nper);
else
    d = s.d(:,ones(1,nper));
    ka = s.ka(:,ones(1,nper));
    kf = s.kf(:,ones(1,nper));
end
eu = real(s.tune);
ea = imag(s.tune);
eu(isnan(eu)) = 0;
ea(isnan(ea)) = 0;
lasta = max([0,find(any(ea ~= 0,1),1,'last')]);
lastu = max([0,find(any(eu ~= 0,1),1,'last')]);
last = max([lasta,lastu]);
if isempty(last)
    return
end
for t = 2 : last
    e = [eu(:,t) + ea(:,t),ea(:,t+1:lasta)];
    k = size(e,2);
    d(:,t) = d(:,t) + H*e(:,1);
    kf(:,t) = kf(:,t) + Rf(:,1:ne*k)*e(:);
    ka(:,t) = ka(:,t) + Ra(:,1:ne*k)*e(:);
end
end
% xxshocktunes().

%**************************************************************************
function P = xxpa2pb(U,P)
% xxpa2pb  Convert MSE matrix for alpha vector to MSE matrix for
% predetermined vector.
Ut = U.';
for i = 1 : size(P,3)
    P(:,:,i) = U*P(:,:,i)*Ut;
    d = diag(P(:,:,i));
    index = d <= 0;
    P(index,:,i) = 0;
    P(:,index,i) = 0;
end
end
% xxpa2pb().

%**************************************************************************
function s = xxomg2sasy(s)
% Convert the structural covariance matrix Omega to reduced-form
% covariance matrices Sa and Sy. Detect `Inf` std deviations and remove
% the corresponding observations.
nb = s.nb;
ny = s.ny;
ne = s.ne;
nper = s.nper;
lastomg = size(s.Omg,3);
tshocks = s.tshocks;
mshocks = s.mshocks;

% Periods where Omg(t) is the same as Omg(t-1).
omgeq = [false,all(s.stdcorr(:,1:end-1) == s.stdcorr(:,2:end),1)];

% Cut off forward expansion.
Ra = s.Ra(:,1:ne);
Ra = Ra(:,tshocks);
Rat = Ra.';

H = s.H(:,mshocks);
Ht = s.H(:,mshocks).';

s.Sa = nan(nb,nb,lastomg);
s.Sy = nan(ny,ny,lastomg);
s.syinf = false(ny,lastomg);

for t = 1 : lastomg
    % If Omg(t) is the same as Omg(t-1), do not compute anything and
    % only copy the previous results.
    if omgeq(t)
        s.Sa(:,:,t) = s.Sa(:,:,t-1);
        s.Sy(:,:,t) = s.Sy(:,:,t-1);
        s.syinf(:,t) = s.syinf(:,t-1);
        continue
    end
    Omg = s.Omg(:,:,t);
    Omga = Omg(tshocks,tshocks);
    s.Sa(:,:,t) = Ra*Omga*Rat;
    Omgy = Omg(mshocks,mshocks);
    omgyinf = isinf(diag(Omgy));
    if ~any(omgyinf)
        % No `Inf` std devs.
        s.Sy(:,:,t) = H*Omgy*Ht;
    else
        % Some std devs are Inf, remove the corresponding observations.
        s.Sy(:,:,t) = ...
            H(:,~omgyinf)*Omgy(~omgyinf,~omgyinf)*Ht(~omgyinf,:);
        s.syinf(:,t) = diag(H(:,omgyinf)*Ht(omgyinf,:)) ~= 0;
    end
end

% Expand `syinf` in 2nd dimension to match the number of periods. This is
% because we use `syinf` to remove observations from `y1` on the whole
% filter range.
if lastomg < nper
    s.syinf(:,end+1:nper) = s.syinf(:,ones(1,nper-lastomg));
end

end
% xxomg2sasy().

%**************************************************************************
function s = xxshocktypes(s,m)
% xxshocktypes  Indices of measurement and transition shocks.
t = m.tzero;
nname = length(m.name);
index = nname*(t-1) + find(m.nametype == 3);
moccur = m.occur(m.eqtntype == 1,index);
toccur = m.occur(m.eqtntype == 2,index);
s.mshocks = any(moccur,1);
s.tshocks = any(toccur,1);
end
% xxshocktypes().