Предположим, что латентный процесс содержит модель 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 и наблюдаемые ответы на filter.

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

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

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

stateIndx = [1 3]; % State indices of interest
FilteredStates = NaN(T,numel(stateIndx));
ForecastedStates = NaN(T,numel(stateIndx));

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

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

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'});