Документация

Предположим, что скрытый процесс включает AR (2) и модель MA (1). Существует 50 периодов и MA (1), процесс выпадает из модели в течение итоговых 25 периодов. Следовательно, уравнение состояния в течение первых 25 периодов

и в течение последних 25 периодов, это

где и являются Гауссовыми со средним значением 0 и стандартным отклонением 1.

Предположение, что ряд запускается в 1,5 и 1, соответственно, генерирует случайную последовательность 50 наблюдений от и.

T = 50;
ARMdl = arima('AR',{0.7,-0.2},'Constant',0,'Variance',1);
MAMdl = arima('MA',0.6,'Constant',0,'Variance',1);
x0 = [1.5 1; 1.5 1];
rng(1);
x = [simulate(ARMdl,T,'Y0',x0(:,1)),...
    [simulate(MAMdl,T/2,'Y0',x0(:,2));nan(T/2,1)]];

Последние 25 значений для моделируемого MA (1) данные являются значениями NaN.

Скрытые процессы измеряются с помощью

в течение первых 25 периодов, и

в течение последних 25 периодов, где является Гауссовым со средним значением 0 и стандартным отклонением 1.

Используйте случайный скрытый процесс состояния (x) и уравнение наблюдения, чтобы сгенерировать наблюдения.

y = 2*nansum(x')'+randn(T,1);

Вместе, скрытые уравнения процесса и наблюдения составляют модель в пространстве состояний. Коэффициенты являются неизвестными параметрами, модель в пространстве состояний

в течение первых 25 периодов,

в течение периода 26, и

в течение последних 24 периодов.

Запишите функцию, которая задает, как параметры в params сопоставляют с матрицами модели в пространстве состояний, значениями начального состояния и типом состояния.

% Copyright 2015 The MathWorks, Inc.

function [A,B,C,D,Mean0,Cov0,StateType] = diffuseAR2MAParamMap(params,T)
%diffuseAR2MAParamMap Time-variant diffuse state-space model parameter
%mapping function
%
% This function maps the vector params to the state-space matrices (A, B,
% C, and D) and the type of state (StateType). From periods 1 to T/2, the
% state model is an AR(2) and an MA(1) model, and the observation model is
% the sum of the two states. From periods T/2 + 1 to T, the state model is
% just the AR(2) model.  The AR(2) model is diffuse.
    A1 = {[params(1) params(2) 0 0; 1 0 0 0; 0 0 0 params(3); 0 0 0 0]};
    B1 = {[1 0; 0 0; 0 1; 0 1]}; 
    C1 = {params(4)*[1 0 1 0]};
    Mean0 = [];
    Cov0 = [];
    StateType = [2 2 0 0];
    A2 = {[params(1) params(2) 0 0; 1 0 0 0]};
    B2 = {[1; 0]};
    A3 = {[params(1) params(2); 1 0]};
    B3 = {[1; 0]}; 
    C3 = {params(5)*[1 0]};
    A = [repmat(A1,T/2,1);A2;repmat(A3,(T-2)/2,1)];
    B = [repmat(B1,T/2,1);B2;repmat(B3,(T-2)/2,1)];
    C = [repmat(C1,T/2,1);repmat(C3,T/2,1)];
    D = 1;
end

Сохраните этот код как файл с именем diffuseAR2MAParamMap на вашем пути MATLAB®.

Создайте рассеянную модель в пространстве состояний путем передачи функционального diffuseAR2MAParamMap как указателя на функцию к dssm.

Mdl = dssm(@(params)diffuseAR2MAParamMap(params,T));

dssm неявно создает рассеянную модель в пространстве состояний. Обычно, вы не можете проверить рассеянные модели в пространстве состояний, которые неявно создаются.

Чтобы оценить параметры, передайте наблюдаемые ответы (y) estimate. Задайте произвольный набор положительных начальных значений для неизвестных параметров.

params0 = 0.1*ones(5,1);
EstMdl = estimate(Mdl,y,params0);

Method: Maximum likelihood (fminunc)
Effective Sample size:             48
Logarithmic  likelihood:     -110.313
Akaike   info criterion:      230.626
Bayesian info criterion:      240.186
      |     Coeff       Std Err   t Stat     Prob  
---------------------------------------------------
 c(1) |  0.44041       0.27687    1.59069  0.11168 
 c(2) |  0.03949       0.29585    0.13349  0.89380 
 c(3) |  0.78364       1.49223    0.52515  0.59948 
 c(4) |  1.64260       0.66737    2.46133  0.01384 
 c(5) |  1.90409       0.49374    3.85648  0.00012 
      |                                            
      |   Final State   Std Dev    t Stat    Prob  
 x(1) | -0.81932       0.46706   -1.75420  0.07940 
 x(2) | -0.29909       0.45939   -0.65107  0.51500 

EstMdl является моделью dssm, содержащей предполагаемые коэффициенты. Поверхности вероятности моделей в пространстве состояний могут содержать локальные максимумы. Поэтому попробуйте несколько начальных значений параметров или рассмотрите использование refine.

Отфильтруйте состояния и получите прогнозы состояния путем передачи EstMdl и наблюдаемых ответов на filter.

[~,~,Output]= filter(EstMdl,y);

Output является T-by-1 массив структур, который содержит отфильтрованные состояния и прогнозы состояния.

Преобразуйте Output в таблицу.

OutputTbl = struct2table(Output);
OutputTbl(1:10,1:5) % Display first ten rows of first five variables

ans =

  10x5 table

    LogLikelihood    FilteredStates    FilteredStatesCov    ForecastedStates    ForecastedStatesCov
    _____________    ______________    _________________    ________________    ___________________

      []              []                 []                   []                   []              
      []              []                 []                   []                   []              
      [-2.3218]       [4x1 double]       [4x4 double]         [4x1 double]         [4x4 double]    
      [-2.4464]       [4x1 double]       [4x4 double]         [4x1 double]         [4x4 double]    
      [-3.8758]       [4x1 double]       [4x4 double]         [4x1 double]         [4x4 double]    
      [-2.5212]       [4x1 double]       [4x4 double]         [4x1 double]         [4x4 double]    
      [-1.9016]       [4x1 double]       [4x4 double]         [4x1 double]         [4x4 double]    
      [-1.9284]       [4x1 double]       [4x4 double]         [4x1 double]         [4x4 double]    
      [-2.4110]       [4x1 double]       [4x4 double]         [4x1 double]         [4x4 double]    
      [-2.6502]       [4x1 double]       [4x4 double]         [4x1 double]         [4x4 double]    

Первые две строки таблицы содержат пустые ячейки или нули, которые соответствуют наблюдениям, требуемым инициализировать рассеянный Фильтр Калмана. Таким образом, SwitchTime равняется 2.

SwitchTime = 2;

Извлеките отфильтрованные и предсказанные состояния из таблицы. Вспомните, что два различных состояния находятся в положениях 1 и 3. Состояния в положениях 2 и 4 помогают задать процессы интереса.

stateIdx = [1 3];                        % State indices of interest

FilteredStates = NaN(T,numel(stateIdx));
ForecastedStates = NaN(T,numel(stateIdx));

for t = (SwitchTime + 1):T
    maxInd = size(Output(t).FilteredStates,1);
    mask = stateIdx <= maxInd;
    FilteredStates(t,mask) = Output(t).FilteredStates(stateIdx(mask),1);
		ForecastedStates(t,mask) = Output(t).ForecastedStates(stateIdx(mask),1);
end

FilteredStates(1:SwitchTime,:) = 0;
ForecastedStates(1:SwitchTime,:) = 0;

Постройте истинные значения состояния, отфильтрованные состояния и прогнозы состояния вместе для каждой модели.

figure
plot(1:T,x(:,1),'-k',1:T,FilteredStates(:,1),':r',...
    1:T,ForecastedStates(:,1),'--g','LineWidth',2);
title('AR(2) State Values')
xlabel('Period')
ylabel('State Value')
legend({'True state values','Filtered state values','State forecasts'});

figure
plot(1:T,x(:,2),'-k',1:T,FilteredStates(:,2),':r',...
    1:T,ForecastedStates(:,2),'--g','LineWidth',2);
title('MA(1) State Values')
xlabel('Period')
ylabel('State Value')
legend({'True state values','Filtered state values','State forecasts'});