Прогнозирование модели пространства состояний с использованием методов Монте-Карло

Открыть Live Script

Этот пример показывает, как предсказать модель пространства состояний с помощью методов Монте-Карло и сравнить прогнозы Монте-Карло с теоретическими прогнозами.

Предположим, что связь между изменением уровня безработицы ( $x_{1, t}$ ) и номинальные темпы роста валового национального продукта (nGNP) ( $x_{3, t}$ ) может быть выражена в следующей форме модели пространства состояний.

$[\begin{array}{cccccccccccccccccccc} x_{1, t} \\ x_{2, t} \\ x_{3, t} \\ x_{4, t} \end{array}] = [\begin{array}{cccccccccccccccccccc} ϕ_{1} & θ_{1} & γ_{1} & 0 \\ 0 & 0 & 0 & 0 \\ γ_{2} & 0 & ϕ_{2} & θ_{2} \\ 0 & 0 & 0 & 0 \end{array}] [\begin{array}{cccccccccccccccccccc} x_{1, t - 1} \\ x_{2, t - 1} \\ x_{3, t - 1} \\ x_{4, t - 1} \end{array}] + [\begin{array}{cccccccccccccccccccc} 1 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & 1 \end{array}] [\begin{array}{cccccccccccccccccccc} u_{1, t} \\ u_{2, t} \end{array}]$

$[\begin{array}{cccccccccccccccccccc} y_{1, t} \\ y_{2, t} \end{array}] = [\begin{array}{cccccccccccccccccccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{array}] [\begin{array}{cccccccccccccccccccc} x_{1, t} \\ x_{2, t} \\ x_{3, t} \\ x_{4, t} \end{array}] + [\begin{array}{cccccccccccccccccccc} σ_{1} & 0 \\ 0 & σ_{2} \end{array}] [\begin{array}{cccccccccccccccccccc} ε_{1, t} \\ ε_{2, t} \end{array}],$

где:

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

Загрузите набор данных Нельсона-Плоссера, который содержит, среди прочего, уровень безработицы и серию 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
y = zeros(T-1,2);                          % Preallocate
y(:,1) = diff(u);
y(:,2) = diff(log(gnpn));

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

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

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

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

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

Задайте модель пространства состояний используя ssm. Проверьте, что спецификация модели соответствует модели пространства состояний.

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

Mdl = 
State-space model type: ssm

State vector length: 4
Observation vector length: 2
State disturbance vector length: 2
Observation innovation vector length: 2
Sample size supported by model: Unlimited
Unknown parameters for estimation: 8

State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...
Unknown parameters: c1, c2,...

State equations:
x1(t) = (c1)x1(t-1) + (c3)x2(t-1) + (c4)x3(t-1) + u1(t)
x2(t) = u1(t)
x3(t) = (c2)x1(t-1) + (c5)x3(t-1) + (c6)x4(t-1) + u2(t)
x4(t) = u2(t)

Observation equations:
y1(t) = x1(t) + (c7)e1(t)
y2(t) = x3(t) + (c8)e2(t)

Initial state distribution:

Initial state means are not specified.
Initial state covariance matrix is not specified.
State types are not specified.

Оцените параметры модели и используйте случайный набор начальных значений параметров для оптимизации. Ограничьте оценку $σ_{1}$ и $σ_{2}$ ко всем положительным, вещественным числам, использующим 'lb' аргумент пары "имя-значение". Для числовой устойчивости задайте Гессиана, когда программное обеспечение вычисляет ковариационную матрицу параметра, используя 'CovMethod' аргумент пары "имя-значение".

rng(1);
params0 = rand(8,1);
[EstMdl,estParams] = estimate(Mdl,isY,params0,...
    'lb',[-Inf -Inf -Inf -Inf -Inf -Inf 0 0],'CovMethod','hessian');

Method: Maximum likelihood (fmincon)
Sample size: 51
Logarithmic  likelihood:      -170.92
Akaike   info criterion:       357.84
Bayesian info criterion:      373.295
      |     Coeff       Std Err    t Stat     Prob  
----------------------------------------------------
 c(1) |  0.06750       0.16548     0.40791  0.68334 
 c(2) | -0.01372       0.05887    -0.23302  0.81575 
 c(3) |  2.71201       0.27039    10.03006   0      
 c(4) |  0.83815       2.84585     0.29452  0.76836 
 c(5) |  0.06273       2.83473     0.02213  0.98235 
 c(6) |  0.05197       2.56875     0.02023  0.98386 
 c(7) |  0.00273       2.40763     0.00113  0.99910 
 c(8) |  0.00016       0.13942     0.00113  0.99910 
      |                                             
      |   Final State   Std Dev     t Stat    Prob  
 x(1) | -0.00000       0.00273    -0.00033  0.99973 
 x(2) |  0.12237       0.92954     0.13164  0.89527 
 x(3) |  0.04049       0.00016   256.76330   0      
 x(4) |  0.01183       0.00016    72.51366   0

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

Фильтрация оценочной модели пространства состояний и извлечение отфильтрованных состояний и их отклонений из последнего периода.

[~,~,Output] = filter(EstMdl,isY);

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

EstMdl1 = EstMdl;
EstMdl1.Mean0 = Output(end).FilteredStates;
EstMdl1.Cov0 = Output(end).FilteredStatesCov;

Моделируйте 5e5 пути наблюдений от установленной модели пространства состояний EstMdl. Задайте, чтобы симулировать наблюдения для каждого периода.

numPaths = 5e5;
SimY = simulate(EstMdl1,10,'NumPaths',numPaths);

SimY является 10-by- 2-by- numPaths массив, содержащий моделируемые наблюдения. Строки SimY соответствуют периодам, столбцы соответствуют наблюдению в модели, а страницы соответствуют путям.

Оцените прогнозные наблюдения и их 95% доверительные интервалы в прогнозном горизонте.

MCFY = mean(SimY,3);
CIFY = quantile(SimY,[0.025 0.975],3);

Оцените теоретические полосы прогноза.

[Y,YMSE] = forecast(EstMdl,10,isY);
Lb = Y - sqrt(YMSE)*1.96;
Ub = Y + sqrt(YMSE)*1.96;

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

figure
h = plot(dates(end-numPeriods-9:end),[isY(end-9:end,1);oosY(:,1)],'-k',...
    dates(end-numPeriods+1:end),MCFY(end-numPeriods+1:end,1),'.-r',...
    dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,1,1),'-b',...
    dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,1,2),'-b',...
    dates(end-numPeriods+1:end),Y(:,1),':c',...
    dates(end-numPeriods+1:end),Lb(:,1),':m',...
    dates(end-numPeriods+1:end),Ub(:,1),':m',...
    'LineWidth',3);
xlabel('Period')
ylabel('Change in the unemployment rate')
legend(h([1,2,4:6]),{'Observations','MC forecasts',...
    '95% forecast intervals','Theoretical forecasts',...
    '95% theoretical intervals'},'Location','Best')
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 7 objects of type line. These objects represent Observations, MC forecasts, 95% forecast intervals, Theoretical forecasts, 95% theoretical intervals.

figure
h = plot(dates(end-numPeriods-9:end),[isY(end-9:end,2);oosY(:,2)],'-k',...
    dates(end-numPeriods+1:end),MCFY(end-numPeriods+1:end,2),'.-r',...
    dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,2,1),'-b',...
    dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,2,2),'-b',...
    dates(end-numPeriods+1:end),Y(:,2),':c',...
    dates(end-numPeriods+1:end),Lb(:,2),':m',...
    dates(end-numPeriods+1:end),Ub(:,2),':m',...
    'LineWidth',3);
xlabel('Period')
ylabel('nGNP growth rate')
legend(h([1,2,4:6]),{'Observations','MC forecasts',...
    '95% MC intervals','Theoretical forecasts','95% theoretical intervals'},...
    'Location','Best')
title('Observed and Forecasted nGNP Growth Rates')

Figure contains an axes. The axes with title Observed and Forecasted nGNP Growth Rates contains 7 objects of type line. These objects represent Observations, MC forecasts, 95% MC intervals, Theoretical forecasts, 95% theoretical intervals.

См. также

estimate | forecast | refine | simulate | ssm

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

Прогнозирование модели пространства состояний с использованием методов Монте-Карло

См. также

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

Подробнее о

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

Поддержка