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