Daniel@0: function [xsmooth, Vsmooth, VVsmooth, loglik] = kalman_smoother(y, A, C, Q, R, init_x, init_V, varargin)
Daniel@0: % Kalman/RTS smoother.
Daniel@0: % [xsmooth, Vsmooth, VVsmooth, loglik] = kalman_smoother(y, A, C, Q, R, init_x, init_V, ...)
Daniel@0: %
Daniel@0: % The inputs are the same as for kalman_filter.
Daniel@0: % The outputs are almost the same, except we condition on y(:, 1:T) (and u(:, 1:T) if specified),
Daniel@0: % instead of on y(:, 1:t).
Daniel@0: 
Daniel@0: [os T] = size(y);
Daniel@0: ss = length(A);
Daniel@0: 
Daniel@0: % set default params
Daniel@0: model = ones(1,T);
Daniel@0: u = [];
Daniel@0: B = [];
Daniel@0: 
Daniel@0: args = varargin;
Daniel@0: nargs = length(args);
Daniel@0: for i=1:2:nargs
Daniel@0:   switch args{i}
Daniel@0:    case 'model', model = args{i+1};
Daniel@0:    case 'u', u = args{i+1};
Daniel@0:    case 'B', B = args{i+1};
Daniel@0:    otherwise, error(['unrecognized argument ' args{i}])
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: xsmooth = zeros(ss, T);
Daniel@0: Vsmooth = zeros(ss, ss, T);
Daniel@0: VVsmooth = zeros(ss, ss, T);
Daniel@0: 
Daniel@0: % Forward pass
Daniel@0: [xfilt, Vfilt, VVfilt, loglik] = kalman_filter(y, A, C, Q, R, init_x, init_V, ...
Daniel@0: 					       'model', model, 'u', u, 'B', B);
Daniel@0: 
Daniel@0: % Backward pass
Daniel@0: xsmooth(:,T) = xfilt(:,T);
Daniel@0: Vsmooth(:,:,T) = Vfilt(:,:,T);
Daniel@0: %VVsmooth(:,:,T) = VVfilt(:,:,T);
Daniel@0: 
Daniel@0: for t=T-1:-1:1
Daniel@0:   m = model(t+1);
Daniel@0:   if isempty(B)
Daniel@0:     [xsmooth(:,t), Vsmooth(:,:,t), VVsmooth(:,:,t+1)] = ...
Daniel@0: 	smooth_update(xsmooth(:,t+1), Vsmooth(:,:,t+1), xfilt(:,t), Vfilt(:,:,t), ...
Daniel@0: 		      Vfilt(:,:,t+1), VVfilt(:,:,t+1), A(:,:,m), Q(:,:,m), [], []);
Daniel@0:   else
Daniel@0:     [xsmooth(:,t), Vsmooth(:,:,t), VVsmooth(:,:,t+1)] = ...
Daniel@0: 	smooth_update(xsmooth(:,t+1), Vsmooth(:,:,t+1), xfilt(:,t), Vfilt(:,:,t), ...
Daniel@0: 		      Vfilt(:,:,t+1), VVfilt(:,:,t+1), A(:,:,m), Q(:,:,m), B(:,:,m), u(:,t+1));
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: VVsmooth(:,:,1) = zeros(ss,ss);
Daniel@0: