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

Прогнозные наблюдения и состояния пяти периодов в будущем. Также получите показатели изменчивости для прогнозов.

numPeriods = 5;
[fY,yMSE,FX,XMSE]= forecast(EstMdl,numPeriods,y);

forecast использует EstMdl.A{end}..., EstMdl.D{end} для прогноза модели пространства состояний. fY и yMSE являются numPeriods-by-1 числовые векторы прогнозируемых наблюдений и отклонений прогнозируемых наблюдений, соответственно. FX и XMSE являются numPeriods-by-2 матрицы прогнозов состояния и отклонений прогнозов состояния. Столбцы указывают состояние, а строки - период. Для всех выходных аргументов последняя строка соответствует последнему прогнозу.

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

figure;
plot(T-10:T,x(T-10:T,1),'-k',T+1:T+numPeriods,FX(:,1),'-r',...
    T-10:T,y(T-10:T),'--g',T+1:T+numPeriods,fY,'--b',...
    T:T+1,[y(T),fY(1);x(T,1),FX(1,1)]',':k','LineWidth',2);
xlabel('Period')
ylabel('States and Observations')
legend({'True state values','State forecasts',...
    'Observed responses','Forecasted responses'});