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

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

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 периодов. В данном примере предположим, что timeseries является стабильным (то есть все параметры являются инвариантными по времени).

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).

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