Класс: ssm
Монте-Карло симуляция государственно-космических моделей
[ моделирует один путь выборки наблюдений (Y,X] = simulate(Mdl,numObs)Y) и состояния (X) из полностью заданной модели состояния-пространства (Mdl). Программное обеспечение моделируется numObs наблюдения и состояния для каждого пути выборки.
[ возвращает смоделированные ответы и состояния с дополнительными опциями, заданными одним или несколькими Y,X] = simulate(Mdl,numObs,Name,Value)Name,Value аргументы пары.
Например, укажите количество путей или значения параметров модели.
Предположим, что скрытый процесс является моделью AR (1). Уравнение состояния
+ ut,
где - гауссов со средним значением 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);
Предположим далее, что скрытый процесс подвержен аддитивной погрешности измерения. Уравнение наблюдения
αt,
где - гауссов со средним значением 0 и стандартным отклонением 0,75. В совокупности скрытый процесс и уравнения наблюдения составляют модель состояния-пространства.
Использовать процесс случайного скрытого состояния (x) и уравнение наблюдения для генерации наблюдений.
y = x + 0.75*randn(T,1);
Укажите четыре матрицы коэффициентов.
A = 0.5; B = 1; C = 1; D = 0.75;
Укажите модель state-space с помощью матриц коэффициентов.
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 моделирует один путь для каждого состояния и наблюдения в модели state-space. Чтобы провести исследование Монте-Карло, укажите, чтобы смоделировать большое количество путей.
Чтобы создать вариации из модели пространства состояний, укажите значения для всех неизвестных параметров.
Явно создайте эту модель состояния-пространства.
где и - независимые гауссовы случайные величины со средним значением 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 значения по столбцам, следующие за порядками A, B, C, D, Mean0 и Cov0. Порядок элементов в params должен соответствовать этому поиску.
Предположим, что взаимосвязь между изменением уровня безработицы (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')
Смоделировать состояния из их совместного условного заднего распределения с учетом ответов с помощью simsmooth.
[1] Дурбин Дж., и С. Дж. Копман. Анализ временных рядов по методам пространства состояний. 2-й ред. Оксфорд: Oxford University Press, 2012.
estimate | filter | forecast | simsmooth | smooth | ssm