simulate

Класс: ssm

Симуляция Монте-Карло моделей пространства состояний

Описание

пример

[Y,X] = simulate(Mdl,numObs) симулирует один выборочный путь наблюдений (Y) и состояний (X) из полностью заданной модели пространства состояний (Mdl). Программа моделирует numObs наблюдения и состояния на путь расчета.

[Y,X] = simulate(Mdl,numObs,Name,Value) возвращает симулированные отклики и состояния с дополнительными опциями, заданными одним или несколькими Name,Value аргументы в виде пар.

Для примера задайте количество путей или значений параметров модели.

пример

[Y,X,U,E] = simulate(___) дополнительно моделируйте нарушения порядка состояния (U) и инновации в области наблюдений (E) с использованием любого из входных параметров в предыдущих синтаксисах.

Входные параметры

расширить все

Стандартная модель пространства состояний, заданная как ssm объект модели, возвращенный ssm или estimate. Стандартная модель пространства состояний имеет конечные начальные элементы ковариационной матрицы состояний. То есть Mdl не может быть dssm объект модели.

Если Mdl не полностью задан (то есть Mdl содержит неизвестные параметры), затем задайте значения для неизвестных параметров используя 'Params' Name,Value аргумент в виде пары. В противном случае программа выдает ошибку.

Количество периодов в пути для генерации вариантов, заданное в виде положительного целого числа.

Если Mdl является изменяющейся во времени моделью, тогда длина вектора камеры, соответствующего матрицам коэффициентов, должна быть по меньшей мере numObs.

Если numObs меньше, чем количество периодов, которые Mdl может поддерживать, тогда программное обеспечение использует только матрицы в первом numObs камеры векторов камеры, соответствующих матрицам коэффициентов.

Типы данных: double

Аргументы в виде пар имя-значение

Задайте необязательные разделенные разделенными запятой парами Name,Value аргументы. Name - имя аргумента и Value - соответствующее значение. Name должны находиться внутри кавычек. Можно задать несколько аргументов в виде пар имен и значений в любом порядке Name1,Value1,...,NameN,ValueN.

Количество путей расчета для генерации вариантов, заданное как разделенная разделенными запятой парами, состоящая из 'NumPaths' и положительное целое число.

Пример: 'NumPaths',1000

Типы данных: double

Значения для неизвестных параметров в модели пространства состояний, заданные как разделенная разделенными запятой парами, состоящая из 'Params' и числовой вектор.

Элементы Params соответствуют неизвестным параметрам в матрицах модели пространства состояний A, B, C, и D, и, опционально, начальное состояние означает Mean0 и ковариационная матрица Cov0.

  • Если вы создали Mdl явным образом (то есть путем задания матриц без функции отображения параметра в матрицу), затем программное обеспечение отображает элементы Params на NaNs в матрицах модели пространства состояний и начальных значениях состояний. Программа ищет NaNs столбцово по порядку A, B, C, D, Mean0, и Cov0.

  • Если вы создали Mdl неявно (то есть путем задания матриц с функцией отображения параметра в матрицу), тогда необходимо задать начальные значения параметров для матриц модели пространства состояний, начальных значений состояний и типов состояний в функции отображения параметра в матрицу.

Если Mdl содержит неизвестные параметры, затем необходимо задать их значения. В противном случае программа игнорирует значение Params.

Типы данных: double

Выходные аргументы

расширить все

Моделируемые наблюдения, возвращенные как матрица или камера числовых векторов.

Если Mdl является инвариантной по времени моделью относительно наблюдений, затем Y является numObs-by- n -by- numPaths массив. То есть каждая строка соответствует периоду, каждый столбец соответствует наблюдению в модели, и каждая страница соответствует пути выборки. Последняя строка соответствует последним моделируемым наблюдениям.

Если Mdl является изменяющейся во времени моделью относительно наблюдений, затем Y является numObs-by- numPaths камера векторов. Y{t,j} содержит вектор длины nt моделируемых наблюдений для t периода j пути выборки. Последняя строка Y содержит последний набор моделируемых наблюдений.

Типы данных: cell | double

Моделируемые состояния, возвращенные как числовая матрица или камера векторов.

Если Mdl является инвариантной по времени моделью относительно состояний, затем X является numObs-by- m -by- numPaths массив. То есть каждая строка соответствует периоду, каждый столбец соответствует состоянию в модели, и каждая страница соответствует пути расчета. Последняя строка соответствует последним моделируемым состояниям.

Если Mdl является изменяющейся во времени моделью относительно состояний, затем X является numObs-by- numPaths камера векторов. X{t,j} содержит вектор длины mt моделируемых состояний для периодического t j пути выборки. Последняя строка X содержит последний набор моделируемых состояний.

Моделируемые нарушения порядка состояния, возвращенные как матрица или матрица камер векторов.

Если Mdl является инвариантной по времени моделью относительно нарушений порядка состояния, затем U является numObs-by- h -by- numPaths массив. То есть каждая строка соответствует периоду, каждый столбец соответствует нарушению порядка состояния в модели, и каждая страница соответствует пути расчета. Последняя строка соответствует последнему моделируемым нарушениям порядка.

Если Mdl является изменяющейся во времени моделью относительно нарушений порядка состояния, затем U является numObs-by- numPaths камера векторов. U{t,j} содержит вектор длины ht моделируемого состояния нарушений порядка для периода t пути выборки j. Последняя строка U содержит последний набор моделируемых нарушений порядка состояния.

Типы данных: cell | double

Моделируемые инновации наблюдений, возвращенные как матрица или камера числовых векторов.

Если Mdl является инвариантной по времени моделью относительно инноваций наблюдений, затем E является numObs-by- h -by- numPaths массив. То есть каждая строка соответствует периоду, каждый столбец соответствует инновациям наблюдений в модели, и каждая страница соответствует пути выборки. Последняя строка соответствует последним моделируемым инновациям в области наблюдений.

Если Mdl является изменяющейся во времени моделью относительно инноваций в области наблюдений, затем E является numObs-by- numPaths камера векторов. E{t,j} содержит вектор длины ht моделируемых инноваций наблюдений для периода t j пути выборки. Последняя строка E содержит последний набор моделируемых наблюдений.

Типы данных: cell | double

Примеры

расширить все

Предположим, что латентный процесс является моделью AR (1). Уравнение состояния

xt=0.5xt-1+ut,

где ut является Гауссовым со средним 0 и стандартным отклонением 1.

Сгенерируйте случайную серию из 100 наблюдений xt, принимая, что серия начинается с 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);

Предположим далее, что скрытый процесс подвержен аддитивной ошибке измерения. Уравнение наблюдения

yt=xt+εt,

где εt является Гауссовым со средним 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'})

Figure contains 2 axes. Axes 1 with title True State Values and Simulated States contains 2 objects of type line. These objects represent True state values, Simulated state values. Axes 2 with title Observed Responses and Simulated responses contains 2 objects of type line. These objects represent Observed responses, Simulated responses.

По умолчанию simulate моделирует один путь для каждого состояния и наблюдения в модели пространства состояний. Чтобы провести исследование Монте-Карло, задайте, чтобы симулировать большое количество путей.

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

Явно создайте эту модель пространства состояний.

xt=ϕxt-1+σ1utyt=xt+σ2εt

где ut и εt являются независимыми Гауссовыми случайными переменными со средним 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';

Figure contains an axes. The axes with title Simulated Responses contains an object of type line.

Программа ищет NaN значения в столбце по порядку A, B, C, D, Mean0 и Cov0. Порядок элементов в params должен соответствовать этому поиску.

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

[x1,tx2,tx3,tx4,t]=[ϕ1θ1γ100000γ20ϕ2θ20000][x1,t-1x2,t-1x3,t-1x4,t-1]+[10100101][u1,tu2,t]

[y1,ty2,t]=[10000010][x1,tx2,tx3,tx4,t]+[σ100σ2][ε1,tε2,t],

где:

  • x1,t - изменение уровня безработицы в момент t.

  • x2,t является фиктивным состоянием для эффекта MA (1) наx1,t.

  • x3,t - скорость роста nGNP в момент t.

  • x4,t является фиктивным состоянием для эффекта MA (1) наx3,t.

  • y1,t - наблюдаемое изменение уровня безработицы.

  • y2,t - наблюдаемая скорость роста nGNP.

  • u1,t и u2,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.

Совет

Симулируйте состояния из их совместного условного апостериорного распределения, учитывая ответы при помощи simsmooth.

Ссылки

[1] Дурбин Дж., и С. Дж. Копман. Анализ временных рядов по методам пространства состояний. 2nd ed. Oxford: Oxford University Press, 2012.