моделировать

Класс: 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.

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

пример

[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 должен появиться в кавычках. Вы можете задать несколько аргументов в виде пар имен и значений в любом порядке, например: Name1, Value1, ..., NameN, ValueN.

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

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

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

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

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

  • Если вы создали Mdl явным образом (то есть, путем определения матриц без функции отображения параметра к матрице), то программное обеспечение сопоставляет элементы Params к NaN s в матрицах модели в пространстве состояний и значениях начального состояния. Программное обеспечение ищет NaN s по столбцам выполняющий приказ 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'})

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

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