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) использование любого из входных параметров в предыдущих синтаксисах.

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

развернуть все

Стандартная модель в пространстве состояний в виде anssm объект модели возвращен 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 соответствуйте неизвестным параметрам в матрицах модели в пространстве состояний ABC, и D, и, опционально, начальное состояние означает Mean0 и ковариационная матрица Cov0.

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

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

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

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

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

развернуть все

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

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

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

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

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

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

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

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

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

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

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

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

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

Если Mdl изменяющаяся во времени модель относительно инноваций наблюдения, затем E numObs- 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'})

По умолчанию, 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';

Программное обеспечение ищет NaN значения, по столбцам выполняющие приказ, 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 серия 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.03005   0      
 c(4) |  0.83816       2.84586     0.29452  0.76836 
 c(5) |  0.06275       2.83467     0.02214  0.98234 
 c(6) |  0.05196       2.56872     0.02023  0.98386 
 c(7) |  0.00274       2.40768     0.00114  0.99909 
 c(8) |  0.00016       0.13942     0.00113  0.99910 
      |                                             
      |   Final State   Std Dev     t Stat    Prob  
 x(1) | -0.00000       0.00274    -0.00033  0.99973 
 x(2) |  0.12237       0.92954     0.13164  0.89527 
 x(3) |  0.04049       0.00016   256.69810   0      
 x(4) |  0.01183       0.00016    72.48544   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')

Советы

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

Ссылки

[1] Дербин Дж. и С. Дж. Купмен. Анализ Временных рядов Методами Пространства состояний. 2-й редактор Оксфорд: Издательство Оксфордского университета, 2012.

Для просмотра документации необходимо авторизоваться на сайте