Прогнозные наблюдения модели пространства состояний, содержащей регрессионый компонент

Открыть Live Script

Этот пример показывает, как оценить регрессионую модель, содержащую регрессионый компонент, а затем предсказать наблюдения из подобранной модели.

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

$\begin{array}{l} [\begin{array}{cccccccccccccccccccc} x_{1, t} \\ x_{2, t} \end{array}] = [\begin{array}{cccccccccccccccccccc} ϕ & θ \\ 0 & 0 \end{array}] [\begin{array}{cccccccccccccccccccc} x_{1, t - 1} \\ x_{2, t - 1} \end{array}] + [\begin{array}{cccccccccccccccccccc} 1 \\ 1 \end{array}] u_{1, t} \\ y_{t} - β Z_{t} = x_{1, t} + σ ε_{t}, \end{array}$

где:

$x_{1, t}$ - изменение уровня безработицы в момент t.
$x_{2, t}$ является фиктивным состоянием для эффекта MA (1).
$y_{1, t}$ наблюдаемое изменение уровня безработицы отклоняется темпами роста nGNP ( $Z_{t}$ ).
$u_{1, t}$ - Гауссов ряд нарушений порядка состояния, имеющих среднее 0 и стандартное отклонение 1.
$ε_{t}$ - гауссов ряд инноваций наблюдений, имеющих среднее 0 и стандартное отклонение $σ$ .

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

load Data_NelsonPlosser

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

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

Хотя этот пример удаляет отсутствующие значения, программное обеспечение может включать серии, содержащие отсутствующие значения в среде фильтра Калмана.

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

numPeriods = 10;                   % Forecast horizon
isY = y(1:end-numPeriods);         % In-sample observations
oosY = y(end-numPeriods+1:end);    % Out-of-sample observations
ISZ = Z(1:end-numPeriods,:);       % In-sample predictors
OOSZ = Z(end-numPeriods+1:end,:);  % Out-of-sample predictors

Задайте матрицы коэффициентов.

A = [NaN NaN; 0 0];
B = [1; 1];
C = [1 0];
D = NaN;

Задайте модель пространства состояний используя ssm.

Mdl = ssm(A,B,C,D);

Оцените параметры модели. Задайте регрессионный компонент и его начальное значение для оптимизации с помощью 'Predictors' и 'Beta0' Аргументы пары "имя-значение", соответственно. Ограничьте оценку $σ$ ко всем положительным, вещественным числам. Для числовой устойчивости задайте Гессиана, когда программное обеспечение вычисляет ковариационную матрицу параметра, используя 'CovMethod' аргумент пары "имя-значение".

params0 = [0.3 0.2 0.1]; % Chosen arbitrarily
[EstMdl,estParams] = estimate(Mdl,isY,params0,'Predictors',ISZ,...
    'Beta0',[0.1 0.2],'lb',[-Inf,-Inf,0,-Inf,-Inf],'CovMethod','hessian');

Method: Maximum likelihood (fmincon)
Sample size: 51
Logarithmic  likelihood:     -87.2409
Akaike   info criterion:      184.482
Bayesian info criterion:      194.141
           |      Coeff       Std Err    t Stat     Prob  
----------------------------------------------------------
 c(1)      |  -0.31780       0.19429    -1.63572  0.10190 
 c(2)      |   1.21242       0.48882     2.48031  0.01313 
 c(3)      |   0.45583       0.63930     0.71301  0.47584 
 y <- z(1) |   1.32407       0.26313     5.03201   0      
 y <- z(2) | -24.48733       1.90115   -12.88024   0      
           |                                              
           |    Final State   Std Dev     t Stat    Prob  
 x(1)      |  -0.38117       0.42842    -0.88971  0.37363 
 x(2)      |   0.23402       0.66222     0.35339  0.72380

EstMdl является ssm модель, и вы можете получить доступ к ее свойствам с помощью записи через точку.

Прогнозные наблюдения за прогнозным горизонтом. EstMdl не хранит набор данных, поэтому необходимо передать его в соответствующих аргументах пары "имя-значение".

[fY,yMSE] = forecast(EstMdl,numPeriods,isY,'Predictors0',ISZ,...
    'PredictorsF',OOSZ,'Beta',estParams(end-1:end));

fY вектор 10 на 1, содержащий прогнозируемые наблюдения и yMSE вектор 10 на 1, содержащий отклонения прогнозируемых наблюдений.

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

ForecastIntervals(:,1) = fY - 1.96*sqrt(yMSE);
ForecastIntervals(:,2) = fY + 1.96*sqrt(yMSE);

figure
h = plot(dates(end-numPeriods-9:end-numPeriods),isY(end-9:end),'-k',...
    dates(end-numPeriods+1:end),oosY,'-k',...
    dates(end-numPeriods+1:end),fY,'--r',...
    dates(end-numPeriods+1:end),ForecastIntervals,':b',...
    dates(end-numPeriods:end-numPeriods+1),...
    [isY(end)*ones(3,1),[oosY(1);ForecastIntervals(1,:)']],':k',...
    'LineWidth',2);
xlabel('Period')
ylabel('Change in the unemployment rate')
legend(h([1,3,4]),{'Observations','Forecasted responses',...
    '95% forecast intervals'})
title('Observed and Forecasted Changes in the Unemployment Rate')

Figure contains an axes. The axes with title Observed and Forecasted Changes in the Unemployment Rate contains 8 objects of type line. These objects represent Observations, Forecasted responses, 95% forecast intervals.

Эта модель, по-видимому, хорошо прогнозирует изменения уровня безработицы.

См. также

estimate | forecast | refine | ssm

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

Прогнозные наблюдения модели пространства состояний, содержащей регрессионый компонент

См. также

Похожие примеры

Подробнее о

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

Поддержка