Discrete time kalman filter with augmented state vector

Question

I am trying to implement the discrete time kalman filter for a state space model with an augmented state vector. That is, the state equation is given by:

x(t) = 0.80x(t-1) + noise

There is only latent process in my example so the x(t) is of dimension 1x1. The measurement equation is given by:

Y(t) = AX(t) + noise

Where capital Y(t) has dimension 3x1. It contains three observed series, Y(t) = [y1(t),y2(t),y3(t)]'. Two of the observed series are occasionally missing in a systematic way:

• y1(t) is observed at all time points.
• y2(t) is observed at every fourth time point.
• y3(t) is observed at every twelfth time point.

The dimension of the measurement matrix A is 3x12 and the dimension of capital X(t) is 12x1. This is a vector stacked with the last twelve realizations of the latent process:

X(t) = [x(t),x(t-1),...,x(t-11)]'.

The reason for capital X(t) being stacked is that y1(t) is related to x(t) and x(t-1). Moreover, y2(t) is related to x(t), x(t-1), x(t-2) and x(t-3). Finally, y3(t) is related to all of the last twelve realizations of the latent process.

The attached matlab code simulates data from this state space model and subsequently runs through a kalman filter with an augmented state space vector, X(t).

My problem is, that the filtered (and the predicted) process differs substantially from the true latent process. This can be seen from the attached figure as well. Clearly, I am doing something wrong but I am unable to see what it is? My kalman filter worked very well, when the state vector was not stacked... I hope someone can help, thanks!

    %-------------------------------------- %
    % AUGMENTED KALMAN FILTER WITH MISSINGS
    % ------------------------------------- %

    clear; close all; clc;

    % 1) SET PARAMETER VALUES %
    % ----------------------- %

    % Number of observations, number of latent processes and observed
    % processes.
    T = 2000;
    NumberOfLatent = 1;
    NumberOfObs = 3;

    % Measurement Matrix, A.
    A = [-0.40, -0.15, zeros(1,10);
          0.25,  0.25, 0.25, 0.25, zeros(1,8);
         ones(1,12)];

    % State transition matrix, Phi. 
    Phi = zeros(12,12);
    Phi(1,1) = 0.80;
    Phi(2:end,1:end-1) = eye(11);

    % Covariance matrix (for the measurement equation).            
    R = diag([2.25; 1.00; 1.00]);

    % State noise and measurement noise.
    MeasNoise = mvnrnd(zeros(1,NumberOfObs),R,T)';
    StateNoise = randn(NumberOfLatent,T+11);

    % One latent process (X) and three observed processes (Y).
    X = zeros(NumberOfObs,T+11);  
    Y = zeros(NumberOfObs,T); 

    % Set initial state.
    X0 = 0;

    % 2) SIMULATE DATA %
    % ---------------- %
    % Before Y (the observed processes) can be simulated we need to simulate 11
    % realisations of X (the latent process).  
    for t = 1:T+11
       if (t == 1)
            X(:,t) = Phi(1,1)*X0 + StateNoise(:,t);
       else
            X(:,t) = Phi(1,1)*X(1,t-1) + StateNoise(:,t);
       end
       if (t>=12)
            Y(:,t-11) = A*X(1,t:-1:t-11)' + MeasNoise(:,t-11);
       end
    end


    % 3) DELETE DATA SUCH THAT THERE ARE MISSINGS %
    % ------------------------------------------- %
    % First row is observed at all time points. The second rows is observed  at
    % every 4th. time point and the third row is observed at every 12th. time
    % point.
    for j = 1:3
        Y(2,j:4:end) = NaN;
    end
    for j = 1:11
        Y(3,j:12:end) = NaN;
    end


    % 4) RUN THE KALMAN FILTER WITH AUGMENTED STATE VECTOR %
    % ---------------------------------------------------- %

    % Initializing matrices.
    AugmentedStates = 12;
    X_predicted = zeros(AugmentedStates,T);                   % Predicted  states.
    X_filtered = zeros(AugmentedStates,T);            % Updated states.
    P_predicted = zeros(AugmentedStates,AugmentedStates,T);   % Predicted variances.
    P_filtered = zeros(AugmentedStates,AugmentedStates,T);   % Filtered variances.
    KG = zeros(AugmentedStates,NumberOfObs,T);       % Kalman gains.
    Q  = eye(AugmentedStates);                       % State noise for  augmented states.
    Res         = zeros(NumberOfObs,T);                       % Residuals.

    % Initial values:
    X_zero = ones(AugmentedStates,1);
    P_zero = eye(AugmentedStates);

    for t = 1:T

        % Initialisation of the Kalman filter.
        if (t == 1)
            X_predicted(:,t) = Phi * X_zero;                          %  Initial predicted state, dimension (AugmentedStates x 1).
            P_predicted(:,:,t) = Phi * P_zero * Phi'+ Q;              % Initial variance, dimension (AugmentedStates x AugmentedStates).
        else
         % #1 Prediction step.
           X_predicted(:,t) = Phi * X_filtered(:,t-1);                    %   Predicted state, dimension (AugmentedStates x 1).
           P_predicted(:,:,t) = Phi * P_filtered(:,:,t-1) * Phi'+ Q;      % Variance of prediction, dimension (AugmentedStates x AugmentedStates).
        end

        % If there is missings the corresponding entries in A and Y is set  to
        % zero.
        if sum(isnan(Y(:,t)))>0
            temp = find(isnan(Y(:,t)));                             
            Y(temp,t) = 0;                                         
            A(temp,:) = 0;              
        end

        % #3 Calculate the Kalman Gains and save them in the matrix KG.
        K = (P_predicted(:,:,t) * A')/(A * P_predicted(:,:,t) * A' + R);         % Kalman Gains, dimension (AugmentedStates x ObservedProcesses).
        KG(:,:,t) = K;  

        % Residuals.
        Res(:,t) = Y(:,t) - A*X_predicted(:,t);                

        % #3 Update step.
        X_filtered(:,t) =  X_predicted(:,t) + K * Res(:,t);                           % Updated state (AugmentedStates x 1).
        P_filtered(:,:,t) = (eye(AugmentedStates) - K * A) * P_predicted(:,:,t);    % Updated variance, dimension (AugmentedStatesstates x AugmentedStates).

         % Measurement Matrix, A, is recreated.
         A = [-0.40, -0.15, zeros(1,10);
               0.25,  0.25, 0.25, 0.25, zeros(1,8);
               ones(1,12)];

     end

Answer 1

Aha. I realized what the bug is. The state transition matrix (Phi) in the matlab code is not correct. That is, in the current matlab code it is:

% State transition matrix, Phi. 
Phi = zeros(12,12);
Phi(1,1) = 0.80;
Phi(2:end,1) = 1

And it should actually look like this:

% State transition matrix, Phi. 
Phi = zeros(12,12);
Phi(1,1) = 0.80;
Phi(2:end,1:end-1) = eye(11)

Without this adjustment the state equation in the filter does not make any sense! I have implemented this in the code. And just to show that it actually works, the plot now looks like this:

... Or at least: The performance of the filter has been improved!