Daniel@0: function [xsmooth, Vsmooth, VVsmooth_future] = smooth_update(xsmooth_future, Vsmooth_future, ...
Daniel@0:     xfilt, Vfilt,  Vfilt_future, VVfilt_future, A, Q, B, u)
Daniel@0: % One step of the backwards RTS smoothing equations.
Daniel@0: % function [xsmooth, Vsmooth, VVsmooth_future] = smooth_update(xsmooth_future, Vsmooth_future, ...
Daniel@0: %    xfilt, Vfilt,  Vfilt_future, VVfilt_future, A, B, u)
Daniel@0: %
Daniel@0: % INPUTS:
Daniel@0: % xsmooth_future = E[X_t+1|T]
Daniel@0: % Vsmooth_future = Cov[X_t+1|T]
Daniel@0: % xfilt = E[X_t|t]
Daniel@0: % Vfilt = Cov[X_t|t]
Daniel@0: % Vfilt_future = Cov[X_t+1|t+1]
Daniel@0: % VVfilt_future = Cov[X_t+1,X_t|t+1]
Daniel@0: % A = system matrix for time t+1
Daniel@0: % Q = system covariance for time t+1
Daniel@0: % B = input matrix for time t+1 (or [] if none)
Daniel@0: % u = input vector for time t+1 (or [] if none)
Daniel@0: %
Daniel@0: % OUTPUTS:
Daniel@0: % xsmooth = E[X_t|T]
Daniel@0: % Vsmooth = Cov[X_t|T]
Daniel@0: % VVsmooth_future = Cov[X_t+1,X_t|T]
Daniel@0: 
Daniel@0: %xpred = E[X(t+1) | t]
Daniel@0: if isempty(B)
Daniel@0:   xpred = A*xfilt;
Daniel@0: else
Daniel@0:   xpred = A*xfilt + B*u;
Daniel@0: end
Daniel@0: Vpred = A*Vfilt*A' + Q; % Vpred = Cov[X(t+1) | t]
Daniel@0: J = Vfilt * A' * inv(Vpred); % smoother gain matrix
Daniel@0: xsmooth = xfilt + J*(xsmooth_future - xpred);
Daniel@0: Vsmooth = Vfilt + J*(Vsmooth_future - Vpred)*J';
Daniel@0: VVsmooth_future = VVfilt_future + (Vsmooth_future - Vfilt_future)*inv(Vfilt_future)*VVfilt_future;
Daniel@0: 
Daniel@0: