Класс: ssm
Симуляция Монте-Карло моделей в пространстве состояний
[Y,X] =
simulate(Mdl,numObs)[Y,X] =
simulate(Mdl,numObs,Name,Value)[Y,X,U,E]
= simulate(___)[ моделирует один демонстрационный путь наблюдений (Y,X] =
simulate(Mdl,numObs)Y) и утверждает (X) от полностью заданной, модели в пространстве состояний (Mdl). Программное обеспечение моделирует наблюдения numObs и состояния на демонстрационный путь.
[ возвращает моделируемые ответы и состояния с дополнительными опциями, заданными одним или несколькими аргументами пары Y,X] =
simulate(Mdl,numObs,Name,Value)Name,Value.
Например, задайте количество значений параметра модели или путей.
Предположим, что скрытый процесс является моделью AR (1). Впоследствии, уравнение состояния
где является Гауссовым со средним значением 0 и стандартным отклонением 1.
Сгенерируйте случайную последовательность 100 наблюдений от , предположение, что ряд запускается в 1,5.
T = 100; ARMdl = arima('AR',0.5,'Constant',0,'Variance',1); x0 = 1.5; rng(1); % For reproducibility x = simulate(ARMdl,T,'Y0',x0);
Предположим далее, что скрытый процесс подвергается аддитивной погрешности измерения. Впоследствии, уравнение наблюдения
где является Гауссовым со средним значением 0 и стандартным отклонением 0.75. Вместе, скрытые уравнения процесса и наблюдения составляют модель в пространстве состояний.
Используйте случайный скрытый процесс состояния (x) и уравнение наблюдения, чтобы сгенерировать наблюдения.
y = x + 0.75*randn(T,1);
Задайте четыре содействующих матрицы.
A = 0.5; B = 1; C = 1; D = 0.75;
Задайте модель в пространстве состояний с помощью содействующих матриц.
Mdl = ssm(A,B,C,D)
Mdl =
State-space model type: ssm
State vector length: 1
Observation vector length: 1
State disturbance vector length: 1
Observation innovation vector length: 1
Sample size supported by model: Unlimited
State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...
State equation:
x1(t) = (0.50)x1(t-1) + u1(t)
Observation equation:
y1(t) = x1(t) + (0.75)e1(t)
Initial state distribution:
Initial state means
x1
0
Initial state covariance matrix
x1
x1 1.33
State types
x1
Stationary
Mdl является моделью ssm. Проверьте, что модель правильно задана с помощью отображения в Командном окне. Программное обеспечение выводит, что процесс состояния является стационарным. Впоследствии, программное обеспечение устанавливает среднее значение начального состояния и ковариацию к среднему значению и отклонению стационарного распределения модели AR (1).
Моделируйте один путь каждое из состояний и наблюдений. Укажите, что пути охватывают 100 периодов.
[simY,simX] = simulate(Mdl,100);
simY 100 1 вектор моделируемых ответов. simX 100 1 вектор моделируемых состояний.
Постройте истинные значения состояния с моделируемыми состояниями. Кроме того, постройте наблюдаемые ответы с моделируемыми ответами.
figure subplot(2,1,1) plot(1:T,x,'-k',1:T,simX,':r','LineWidth',2) title({'True State Values and Simulated States'}) xlabel('Period') ylabel('State') legend({'True state values','Simulated state values'}) subplot(2,1,2) plot(1:T,y,'-k',1:T,simY,':r','LineWidth',2) title({'Observed Responses and Simulated responses'}) xlabel('Period') ylabel('Response') legend({'Observed responses','Simulated responses'})
По умолчанию simulate моделирует один путь для каждого состояния и наблюдения в модели в пространстве состояний. Чтобы провести исследование Монте-Карло, задайте, чтобы моделировать большое количество путей.
Чтобы сгенерировать варьируемые величины от модели в пространстве состояний, задайте значения для всех неизвестных параметров.
Явным образом создайте эту модель в пространстве состояний.
где и независимые Гауссовы случайные переменные со средним значением 0 и отклонением 1. Предположим, что среднее значение начального состояния и отклонение равняются 1, и что состояние является стационарным процессом.
A = NaN; B = NaN; C = 1; D = NaN; mean0 = 1; cov0 = 1; stateType = 0; Mdl = ssm(A,B,C,D,'Mean0',mean0,'Cov0',cov0,'StateType',stateType);
Моделируйте 100 ответов от Mdl. Укажите, что авторегрессивный коэффициент 0.75, стандартное отклонение воздействия состояния 0.5, и инновационное стандартное отклонение наблюдения 0.25.
params = [0.75 0.5 0.25]; y = simulate(Mdl,100,'Params',params); figure; plot(y); title 'Simulated Responses'; xlabel 'Period';
Программное обеспечение ищет значения NaN, по столбцам выполняющие приказ, B, C, D, Mean0 и Cov0. Порядок элементов в params должен соответствовать этому поиску.
Предположим что отношение между изменением в уровне безработицы () и темп роста номинального валового национального продукта (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.83817 2.84587 0.29452 0.76836
c(5) | 0.06274 2.83470 0.02213 0.98234
c(6) | 0.05196 2.56873 0.02023 0.98386
c(7) | 0.00272 2.40773 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.73046 0
x(4) | 0.01183 0.00016 72.51551 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 является 2 10 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')
Моделируйте состояния от их объединенного условного апостериорного распределения, учитывая ответы при помощи simsmooth.
[1] Дербин Дж. и С. Дж. Купмен. Анализ Временных рядов Методами Пространства состояний. 2-й редактор Оксфорд: Издательство Оксфордского университета, 2012.
estimate | filter | forecast | simsmooth | smooth | ssm