В этом примере показано, как предсказать модель в пространстве состояний с помощью методов Монте-Карло, и сравнить прогнозы Монте-Карло с теоретическими прогнозами.
Предположим что отношение между изменением в уровне безработицы () и темп роста номинального валового национального продукта (nGNP) () может быть описан в следующем, форме модели в пространстве состояний.
где:
изменение в уровне безработицы во время t.
фиктивное состояние для MA (1) эффект на .
nGNP темп роста во время t.
фиктивное состояние для MA (1) эффект на .
наблюдаемое изменение в уровне безработицы.
наблюдаемый nGNP темп роста.
и серия Gaussian воздействий состояния, имеющих среднее значение 0 и стандартное отклонение 1.
серия Gaussian инноваций наблюдения, имеющих среднее значение 0 и стандартное отклонение .
серия Gaussian инноваций наблюдения, имеющих среднее значение 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 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.
Оцените параметры модели и используйте случайный набор начальных значений параметров для оптимизации. Ограничьте оценку и ко всем положительным, вещественным числам с помощью '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
- 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 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')
estimate
| forecast
| refine
| simulate
| ssm