Предположим, что скрытый процесс включает 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*sum(x','omitnan')'+randn(T,1);

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

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

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

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

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

% Copyright 2015 The MathWorks, Inc.

function [A,B,C,D,Mean0,Cov0,StateType] = AR2MAParamMap(params,T)
%AR2MAParamMap Time-variant state-space model parameter mapping function
%
% This function maps the vector params to the state-space matrices (A, B,
% C, and D), the initial state value and the initial state variance (Mean0
% and Cov0), 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.
    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 = ones(4,1);
    Cov0 = 10*eye(4);
    StateType = [0 0 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

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

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

Mdl = ssm(@(params)AR2MAParamMap(params,T));

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

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

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

Method: Maximum likelihood (fminunc)
Sample size: 50
Logarithmic  likelihood:     -114.957
Akaike   info criterion:      239.913
Bayesian info criterion:      249.473
      |     Coeff       Std Err   t Stat     Prob  
---------------------------------------------------
 c(1) |  0.47870       0.26634    1.79733  0.07229 
 c(2) |  0.00809       0.27179    0.02975  0.97626 
 c(3) |  0.55735       0.80958    0.68844  0.49118 
 c(4) |  1.62679       0.41622    3.90848  0.00009 
 c(5) |  1.90021       0.49563    3.83391  0.00013 
      |                                            
      |   Final State   Std Dev    t Stat    Prob  
 x(1) | -0.81229       0.46815   -1.73511  0.08272 
 x(2) | -0.31449       0.45918   -0.68490  0.49341 

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

Сглаживайте состояния и оцените ковариационные матрицы отклонения сглаживавших состояний путем передачи EstMdl и наблюдаемые ответы на smooth.

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

Output T- 1 массив структур, содержащий сглаживавшие состояния и их ковариационные матрицы отклонения, среди прочего.

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

stateIndx = [1 3]; % State Indices of interest
SmoothedStates = NaN(T,numel(stateIndx));
SmoothedStatesCov = NaN(T,numel(stateIndx));

for t = 1:T
    maxInd1 = size(Output(t).SmoothedStates,1);
		maxInd2 = size(Output(t).SmoothedStatesCov,1);
    mask1 = stateIndx <= maxInd1;
		mask2 = stateIndx <= maxInd2;
    SmoothedStates(t,mask1) = ...
        Output(t).SmoothedStates(stateIndx(mask1),1);
		SmoothedStatesCov(t,mask2) = ...
        diag(Output(t).SmoothedStatesCov(stateIndx(mask2),...
        stateIndx(mask2)));
end

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

AR2SSCIlb = SmoothedStates(:,1) - 1.95*sqrt(SmoothedStatesCov(:,1));
AR2SSCIub = SmoothedStates(:,1) + 1.95*sqrt(SmoothedStatesCov(:,1));
AR2SSIntervals = [AR2SSCIlb AR2SSCIub];

MA1SSCIlb = SmoothedStates(:,2) - 1.95*sqrt(SmoothedStatesCov(:,2));
MA1SSCIub = SmoothedStates(:,2) + 1.95*sqrt(SmoothedStatesCov(:,2));
MA1SSIntervals = [MA1SSCIlb MA1SSCIub];

figure
plot(1:T,x(:,1),'-k',1:T,SmoothedStates(:,1),':r',...
		1:T,AR2SSIntervals,'--b','LineWidth',2);
title('AR(2) State Values')
xlabel('Period')
ylabel('State Value')
legend({'True state values','Smoothed state values',...
    '95% Confidence Intervals'});

figure
plot(1:T,x(:,2),'-k',1:T,SmoothedStates(:,2),':r',...
		1:T,MA1SSIntervals,'--b','LineWidth',2);
title('MA(1) State Values')
xlabel('Period')
ylabel('State Value')
legend({'True state values','Smoothed state values',...
    '95% Confidence Intervals'});