Предскажите модель в пространстве состояний Используя методы Монте-Карло

Скрипт Open 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}$ серия Gaussian воздействий состояния, имеющих среднее значение 0 и стандартное отклонение 1.
$ε_{1, t}$ серия Gaussian инноваций наблюдения, имеющих среднее значение 0 и стандартное отклонение $σ_{1}$ .
$ε_{2, t}$ серия Gaussian инноваций наблюдения, имеющих среднее значение 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.83816       2.84586     0.29452  0.76836 
 c(5) |  0.06274       2.83469     0.02213  0.98234 
 c(6) |  0.05196       2.56873     0.02023  0.98386 
 c(7) |  0.00272       2.40766     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.00272    -0.00033  0.99973 
 x(2) |  0.12237       0.92954     0.13164  0.89527 
 x(3) |  0.04049       0.00016   256.70438   0      
 x(4) |  0.01183       0.00016    72.50104   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- 2- 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 object. The axes object 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 object. The axes object 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.

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

Предскажите модель в пространстве состояний Используя методы Монте-Карло

Смотрите также

Связанные примеры

Больше о

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

Поддержка