带增广状态向量的离散时间卡尔曼滤波

Question

我试图实现一个状态空间模型的离散时间卡尔曼滤波，一个增广的状态向量。也就是说，状态方程由以下几个方面提供：

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

在我的例子中只有潜在过程，所以x(t)是维数1x1的。计量方程如下：

Y(t) = AX(t) + noise

其中capital Y(t)具有维数3x1。它包含三个观测序列：Y(t) = y1(t)、y2(t)、y3(t)‘。观察到的两个系列偶尔会有系统地缺失：

• y1(t)在所有时间点都被观测到。
• y2(t)每四个时间点观察一次。
• y3(t)每12个时间点观察一次。

度量矩阵A的维数为3x12，资本X(t)的维数为12x1。这是一个向量，与最后12个潜在过程的实现叠加在一起：

X(t) = x(t)，x(t-1)，.，x(t-11)'.

资本市场( capital X(t) )之所以被叠加，是因为y1(t)与x(t)和x(t-1)有关。此外，y2(t)与x(t)、x(t-1)、x(t-2)和xE 267(t-3)有关。最后，y3(t)与最后12个潜在过程的实现有关。

所附的matlab代码模拟了该状态空间模型的数据，并随后通过带有增广状态空间矢量的卡尔曼滤波器( X(t) )运行。

我的问题是，即过滤(和预测)过程与真正的潜在过程有很大的不同。从附图中也可以看出这一点。很明显，我做错了什么，但我看不出来是什么？当状态向量不叠加时，我的卡尔曼滤波器工作得很好.我希望有人能帮忙，谢谢！

    %-------------------------------------- %
    % 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

啊哈。我意识到窃听器是什么。matlab代码中的状态转移矩阵(Phi)是不正确的。也就是说，在当前的matlab代码中，它是：

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

实际上应该是这样的：

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

如果没有这种调整，滤波器中的状态方程就没有任何意义了！我已经在代码中实现了这一点。为了证明它真的有效，这个情节现在看起来是这样的：

​

​

..。或者至少:过滤器的性能已经提高了！