Этот пример показывает, как сгенерировать данные из известной модели, подгонять модель рассеянного пространства состояний к данным, а затем сглаживать состояния.
Предположим, что латентный процесс содержит модель AR (2) и MA (1). Существует 50 периодов, и процесс MA (1) выпадает из модели на последние 25 периодов. Следовательно, уравнение состояния для первых 25 периодов является
и за последние 25 периодов это
где
и
являются Гауссовыми со средним 0 и стандартным отклонением 1.
Принимая, что серия начинается с 1,5 и 1, соответственно, генерируют случайную серию из 50 наблюдений и
T = 50; ARMdl = arima('AR',{0.7,-0.2},'Constant',0,'Variance',1); MAMdl = arima('MA',0.6,'Constant',0,'Variance',1); x0 = [1.5 1; 1.5 1]; rng(1); x = [simulate(ARMdl,T,'Y0',x0(:,1)),... [simulate(MAMdl,T/2,'Y0',x0(:,2));nan(T/2,1)]];
Последние 25 значений для моделируемых данных MA (1) NaN значения.
Латентные процессы измеряются с помощью
за первые 25 периодов, и
за последние 25 периодов, где
Гауссов со средним 0 и стандартным отклонением 1.
Используйте процесс случайного скрытого состояния (x) и уравнение наблюдения для генерации наблюдений.
y = 2*sum(x','omitnan')'+randn(T,1);
Вместе латентный процесс и уравнения наблюдений составляют модель пространства состояний. Если коэффициенты являются неизвестными параметрами, модель пространства состояний
![$$\begin{array}{l}
\left[ {\begin{array}{*{20}{c}}
{{x_{1,t}}}\\
{{x_{2,t}}}\\
{{x_{3,t}}}\\
{{x_{4,t}}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{\phi _1}}&{{\phi _2}}&0&0\\
1&0&0&0\\
0&0&0&{{\theta _1}}\\
0&0&0&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_{1,t - 1}}}\\
{{x_{2,t - 1}}}\\
{{x_{3,t - 1}}}\\
{{x_{4,t - 1}}}
\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
1&0\\
0&0\\
0&1\\
0&1
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{u_{1,t}}}\\
{{u_{2,t}}}
\end{array}} \right]\\
{y_t} = a({x_{1,t}} + {x_{3,t}}) + {\varepsilon _t}
\end{array}$$](../examples/econ/win64/SmoothTimeVaryingDiffuseStateSpaceModelExample_eq06384238688281611023.png)
за первые 25 периодов,
![$$\begin{array}{l}
\left[ {\begin{array}{*{20}{c}}
{{x_{1,t}}}\\
{{x_{2,t}}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{\phi _1}}&{{\phi _2}}&0&0\\
1&0&0&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_{1,t - 1}}}\\
{{x_{2,t - 1}}}\\
{{x_{3,t - 1}}}\\
{{x_{4,t - 1}}}
\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
1\\
0
\end{array}} \right]{u_{1,t}}\\
{y_t} = b{x_{1,t}} + {\varepsilon _t}
\end{array}$$](../examples/econ/win64/SmoothTimeVaryingDiffuseStateSpaceModelExample_eq06140899585903970043.png)
для периода 26, и
![$$\begin{array}{l}
\left[ {\begin{array}{*{20}{c}}
{{x_{1,t}}}\\
{{x_{2,t}}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{\phi _1}}&{{\phi _2}}\\
1&0
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{x_{1,t - 1}}}\\
{{x_{2,t - 1}}}
\end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
1\\
0
\end{array}} \right]{u_{1,t}}\\
{y_t} = b{x_{1,t}} + {\varepsilon _t}
\end{array}$$](../examples/econ/win64/SmoothTimeVaryingDiffuseStateSpaceModelExample_eq09337055536091307358.png)
за последние 24 периода.
Написание функции, которая задает, как параметры в params сопоставить с матрицами модели пространства состояний, начальными значениями состояний и типом состояния.
% Copyright 2015 The MathWorks, Inc. function [A,B,C,D,Mean0,Cov0,StateType] = diffuseAR2MAParamMap(params,T) %diffuseAR2MAParamMap Time-variant diffuse state-space model parameter %mapping function % % This function maps the vector params to the state-space matrices (A, B, % C, and D) and the type of state (StateType). From periods 1 to T/2, the % state model is an AR(2) and an MA(1) model, and the observation model is % the sum of the two states. From periods T/2 + 1 to T, the state model is % just the AR(2) model. The AR(2) model is diffuse. A1 = {[params(1) params(2) 0 0; 1 0 0 0; 0 0 0 params(3); 0 0 0 0]}; B1 = {[1 0; 0 0; 0 1; 0 1]}; C1 = {params(4)*[1 0 1 0]}; Mean0 = []; Cov0 = []; StateType = [2 2 0 0]; A2 = {[params(1) params(2) 0 0; 1 0 0 0]}; B2 = {[1; 0]}; A3 = {[params(1) params(2); 1 0]}; B3 = {[1; 0]}; C3 = {params(5)*[1 0]}; A = [repmat(A1,T/2,1);A2;repmat(A3,(T-2)/2,1)]; B = [repmat(B1,T/2,1);B2;repmat(B3,(T-2)/2,1)]; C = [repmat(C1,T/2,1);repmat(C3,T/2,1)]; D = 1; end
Сохраните этот код как файл с именем diffuseAR2MAParamMap на пути MATLAB ®.
Создайте модель рассеянного пространства состояний путем передачи функции diffuseAR2MAParamMap как указатель на функцию, чтобы dssm.
Mdl = dssm(@(params)diffuseAR2MAParamMap(params,T));
dssm неявно создает модель рассеянного пространства состояний. Обычно вы не можете проверить модели рассеянного пространства состояний, которые неявно созданы.
Чтобы оценить параметры, передайте наблюдаемые отклики (y) к estimate. Задайте произвольный набор положительных начальных значений для неизвестных параметров.
params0 = 0.1*ones(5,1); EstMdl = estimate(Mdl,y,params0);
Method: Maximum likelihood (fminunc)
Effective Sample size: 48
Logarithmic likelihood: -110.313
Akaike info criterion: 230.626
Bayesian info criterion: 240.186
| Coeff Std Err t Stat Prob
---------------------------------------------------
c(1) | 0.44041 0.27687 1.59069 0.11168
c(2) | 0.03949 0.29585 0.13349 0.89380
c(3) | 0.78364 1.49222 0.52515 0.59948
c(4) | 1.64260 0.66736 2.46134 0.01384
c(5) | 1.90409 0.49374 3.85648 0.00012
|
| Final State Std Dev t Stat Prob
x(1) | -0.81932 0.46706 -1.75420 0.07940
x(2) | -0.29909 0.45939 -0.65107 0.51500
EstMdl является dssm модель, содержащая оцененные коэффициенты. Поверхности вероятностей моделей пространства состояний могут содержать локальные максимумы. Поэтому попробуйте несколько начальных значений параметров или рассмотрите использование refine.
Сглаживайте состояния и получайте сглаженную ковариационную матрицу состояния за период путем прохождения EstMdl и наблюдаемые ответы на smooth.
[~,~,Output]= smooth(EstMdl,y);
Output является T-by-1 массив структур, который содержит сглаженные состояния.
Преобразование Output в таблицу.
OutputTbl = struct2table(Output);
OutputTbl(1:10,1:4) % Display first ten rows of first four variables
ans =
10x4 table
LogLikelihood SmoothedStates SmoothedStatesCov SmoothedStateDisturb
_____________ ______________ _________________ ____________________
{0x0 double} {0x0 double} {0x0 double} {0x0 double}
{0x0 double} {0x0 double} {0x0 double} {0x0 double}
{[ -2.3218]} {4x1 double} {4x4 double} {2x1 double}
{[ -2.4464]} {4x1 double} {4x4 double} {2x1 double}
{[ -3.8758]} {4x1 double} {4x4 double} {2x1 double}
{[ -2.5212]} {4x1 double} {4x4 double} {2x1 double}
{[ -1.9016]} {4x1 double} {4x4 double} {2x1 double}
{[ -1.9284]} {4x1 double} {4x4 double} {2x1 double}
{[ -2.4110]} {4x1 double} {4x4 double} {2x1 double}
{[ -2.6502]} {4x1 double} {4x4 double} {2x1 double}
Первые две строки таблицы содержат пустые камеры или нули, которые соответствуют наблюдениям, необходимым для инициализации диффузного фильтра Калмана. То есть SwitchTime равен 2.
SwitchTime = 2;
Извлеките сглаженные состояния из Output, и вычислите их 95% индивидуальных доверительных интервалов типа Уолда. Напомним, что два разных состояния находятся в позициях 1 и 3. Состояния в позициях 2 и 4 помогают определить интересующие процессы.
stateIdx = [1 3]; % State indices of interest SmoothedStates = nan(T,numel(stateIdx)); CI = nan(T,2,numel(stateIdx)); for t = (SwitchTime + 1):T maxInd = size(Output(t).SmoothedStates,1); mask = stateIdx <= maxInd; SmoothedStates(t,mask) = Output(t).SmoothedStates(stateIdx(mask),1); CovX = Output(t).SmoothedStatesCov(stateIdx(mask),stateIdx(mask)); CI(t,:,1) = SmoothedStates(t,1) + 1.96*sqrt(CovX(1,1))*[-1 1]; if (max(stateIdx(mask)) > 1) CI(t,:,2) = SmoothedStates(t,2) + 1.96*sqrt(CovX(2,2))*[-1 1]; end end SmoothedStates(1:SwitchTime,:) = 0; CI(1:SwitchTime,:,:) = 0;
Постройте график значений истинного состояния, сглаженных состояний и 95% доверительных интервалов сглаженного состояния для каждой модели.
figure plot(1:T,x(:,1),'b',1:T,SmoothedStates(:,1),'k',1:T,CI(:,:,1),'--r'); title('AR(2) State Values') xlabel('Period') ylabel('State Value') legend({'True state values','Smoothed state values','95% CI'}); figure plot(1:T,x(:,2),'b',1:T,SmoothedStates(:,2),'k',1:T,CI(:,:,2),'--r'); title('MA(1) State Values') xlabel('Period') ylabel('State Value') legend({'True state values','Smoothed state values','95% CI'});
