Сглаживайте изменяющуюся во времени модель рассеянного пространства состояний

Открыть скрипт

Этот пример показывает, как сгенерировать данные из известной модели, подгонять модель рассеянного пространства состояний к данным, а затем сглаживать состояния.

Предположим, что латентный процесс содержит модель AR (2) и MA (1). Существует 50 периодов, и процесс MA (1) выпадает из модели на последние 25 периодов. Следовательно, уравнение состояния для первых 25 периодов является

$$\begin{array}{l}
{x_{1,t}} = 0.7{x_{1,t - 1}} - 0.2{x_{1,t - 2}} + {u_{1,t}}\\
{x_{2,t}} = {u_{2,t}} + 0.6{u_{2,t - 1}},
\end{array}$$

и за последние 25 периодов это

$${x_{1,t}} = 0.7{x_{1,t - 1}} - 0.2{x_{1,t - 2}} + {u_{1,t}},$$

где$u_{1,t}$ и$u_{2,t}$ являются Гауссовыми со средним 0 и стандартным отклонением 1.

Принимая, что серия начинается с 1,5 и 1, соответственно, генерируют случайную серию из 50 наблюдений и$x_{1,t}$ $x_{2,t}.$

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 значения.

Латентные процессы измеряются с помощью

$${y_t} = 2\left( {{x_{1,t}} + {x_{2,t}}} \right) + {\varepsilon _t},$$

за первые 25 периодов, и

$${y_t} = 2{x_{1,t}} + {\varepsilon _t}$$

за последние 25 периодов, где$\varepsilon_t$ Гауссов со средним 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}$$

за первые 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}$$

для периода 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}$$

за последние 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'});

См. также

| |

Похожие примеры

Подробнее о

Документация по Econometrics Toolbox

Поддержка

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