Предположим, что интерес представляет линейная зависимость между изменением наблюдаемого уровня безработицы и номинальным темпом роста валового национального продукта (нВНП). Предположим далее, что вы хотите выбрать между моделью AR (1) или AR (2) для первой разницы уровня безработицы (т.е. штата). То есть

Объявление функций rwParamMap.m и rwAR2ParamMap.m, которые определяют, каким образом параметры в params сопоставить с матрицами модели state-space, начальными значениями состояния и типом состояния. Сохраните их в рабочей папке.

function [A,B,C,D,Mean0,Cov0,StateType,deflateY] = rwParamMap(params,y,Z)
%rwParamMap Parameter-to-matrix mapping function for rolling window example
%using ssm
%   The state space model specified by rwParamMap contains a stationary
%   AR(1) state, the observation model includes a regression component, and
%   the variances of the innovation and disturbances are 1. The response y
%   is deflated by the regression component specified by the predictor
%   variables x.
A = params(1);
B = 1;
C = 1;
D = 1;
Mean0 = [];
Cov0 = [];
StateType = 0;
deflateY = y - params(2)*Z;
end

function [A,B,C,D,Mean0,Cov0,StateType,deflateY] = rwAR2ParamMap(params,y,Z)
%rwParamMap Parameter-to-matrix mapping function for rolling window example
%using ssm and specifying an ARMA state model
%   The state space model specified by rwParamMap contains a stationary
%   AR(2) state, the observation model includes a regression component, and
%   the variances of the innovation and disturbances are 1. The response y
%   is deflated by the regression component specified by the predictor
%   variables x.
A = [params(1) params(2); 1 0];
B = [1; 0];
C = [1 0];
D = 1;
Mean0 = [];
Cov0 = [];
StateType = [0 0];
deflateY = y - params(3)*Z;
end

Загрузите набор данных Нельсона-Плоссера, который содержит, среди прочего, уровень безработицы и ряды nGNP.

load Data_NelsonPlosser

Выполните предварительную обработку данных, взяв натуральный логарифм ряда nGNP и первую разность каждого. Также снимите пусковую NaN значения из каждой серии.

isNaN = any(ismissing(DataTable),2); % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
ur = DataTable.UR(~isNaN);           % Sample size
Z = diff(log(gnpn));
y = diff(ur);
T = size(y,1);

При неявном определении модели состояния-пространства и данных ответа и предиктора (т. е. y и Z) существуют в MATLAB ® Workspace, то программное обеспечение создает связь из функции отображения параметров в матрицу этих рядов. Если данные не существуют в MATLAB ® Workspace, программа создает модель, но при оценке модели необходимо предоставить данные с использованием соответствующих аргументов пары имя-значение.

Поэтому для выполнения анализа скользящего окна при неявном определении модели «состояние-пространство» и наличии регрессионного компонента необходимо указать модель «состояние-пространство», указывающую индексы данных, подлежащих анализу для каждого окна. Выполните анализ скользящего окна моделируемых данных. Пусть длина перекатываемого окна (mбыть 40 периодами и прогнозным горизонтом (h) 10 периодов. Для этого примера предположим, что временные ряды стабильны (т.е. все параметры инвариантны по времени).

m = 40;
N = T - m + 1;      % Number of rolling windows
h = 10;

fError1 = nan(N,h); % Preallocation
fError2 = nan(N,h);

for j = 1:N;
    idxRW = j:(m + j - h - 1);
		idxFH = (m + j - h):(m + j - 1);
    Mdl1 = ssm(@(c)rwParamMap(c,y(idxRW),Z(idxRW)));
    Mdl2 = ssm(@(c)rwAR2ParamMap(c,y(idxRW),Z(idxRW)));
    [EstMdl1,estParams1] = estimate(Mdl1,y(idxRW),[0.5 -20]',...
        'Display','off');
    [EstMdl2,estParams2] = estimate(Mdl2,y(idxRW),[0.5 0.1 -20]',...
       'Display','off');
    fY1 = forecast(EstMdl1,h,y(idxRW),'Predictors0',Z(idxRW),...
        'PredictorsF',Z(idxFH),'Beta',estParams1(end));
    fY2 = forecast(EstMdl2,h,y(idxRW),'Predictors0',Z(idxRW),...
        'PredictorsF',Z(idxFH),'Beta',estParams2(end));
		fError1(j,:) = y(idxFH) - fY1;
		fError2(j,:) = y(idxFH) - fY2;
end

Вычислите RMSE для каждого перспективного прогноза и сравните их для каждой модели.

fRMSE1 = sqrt(mean(fError1.^2));
fRMSE2 = sqrt(mean(fError2.^2));
fRMSE1 < fRMSE2

ans =

  1x10 logical array

   1   1   1   0   0   0   0   0   0   0

В целом, прогностическая эффективность модели AR (2) лучше, чем модели AR (1).

В качестве альтернативы можно сравнить прогностическую производительность моделей с помощью теста Дибольда-Мариано. Для получения дополнительной информации см. [1].