Документация

Предположим, что линейное соотношение между изменением в наблюдаемом уровне безработицы и темпом роста номинального валового национального продукта (nGNP) представляет интерес. Предположим далее, что вы хотите выбрать между 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).

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