Предположим, что скрытый процесс является случайным обходом. Впоследствии, уравнение состояния

где $$u_t$$ является Гауссовым со средним значением 0 и стандартным отклонением 1.

Сгенерируйте случайную последовательность 100 наблюдений от $$x_t$$, предположение, что ряд запускается в 1,5.

T = 100;
x0 = 1.5;
rng(1); % For reproducibility
u = randn(T,1);
x = cumsum([x0;u]);
x = x(2:end);

Предположим далее, что скрытый процесс подвергается аддитивной погрешности измерения. Впоследствии, уравнение наблюдения

где $${\epsilon }_t$$ является Гауссовым со средним значением 0 и стандартным отклонением 0.75. Вместе, скрытые уравнения процесса и наблюдения составляют модель в пространстве состояний.

Используйте случайный скрытый процесс состояния (x) и уравнение наблюдения, чтобы сгенерировать наблюдения.

y = x + 0.75*randn(T,1);

Задайте четыре содействующих матрицы.

A = 1;
B = 1;
C = 1;
D = 0.75;

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

Mdl = dssm(A,B,C,D,'StateType',2)

Mdl = 
State-space model type: dssm

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) = 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  Inf 

State types
    x1   
 Diffuse 

Mdl является моделью dssm. Проверьте, что модель правильно задана с помощью отображения в Командном окне.

Сглаженные состояния в течение периодов 1 - 100. Постройте истинные значения состояния и сглаживавшие оценки состояния.

SmoothedX = smooth(Mdl,y);

figure
plot(1:T,x,'-k',1:T,SmoothedX,':r','LineWidth',2)
title({'State Values'})
xlabel('Period')
ylabel('State')
legend({'True state values','Smoothed state values'})

Истинные значения и сглаживавшие оценки являются приблизительно тем же самым.

Предположим, что линейное соотношение между уровнем безработицы и номинальным валовым национальным продуктом (nGNP) представляет интерес. Предположим далее, что уровень безработицы является серией AR (1). Символически, и в форме пространства состояний, модель

где:

Загрузите набор данных Нельсона-Плоссера, который содержит уровень безработицы и nGNP ряд, среди прочего.

load Data_NelsonPlosser

Предварительно обработайте данные путем взятия натурального логарифма nGNP ряда и удаления стартовых значений NaN из каждого ряда.

isNaN = any(ismissing(DataTable),2);       % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
y = diff(DataTable.UR(~isNaN));
T = size(gnpn,1);                          % The sample size
Z = price2ret(gnpn);

Этот пример продолжает использовать ряд без значений NaN. Однако с помощью среды Фильтра Калмана, программное обеспечение может разместить ряд, содержащий отсутствующие значения.

Задайте содействующие матрицы.

A = NaN;
B = NaN;
C = 1;

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

StateType = 2;
Mdl = dssm(A,B,C,'StateType',StateType);

Оцените параметры. Задайте компонент регрессии и его начальное значение для оптимизации с помощью аргументов пары "имя-значение" 'Predictors' и 'Beta0', соответственно. Отобразите оценки и всю информацию о диагностике оптимизации. Ограничьте оценку $$\sigma$$ ко всем положительным, вещественным числам.

params0 = [0.3 0.2]; % Initial values chosen arbitrarily
Beta0 = 0.1;
[EstMdl,estParams] = estimate(Mdl,y,params0,'Predictors',Z,'Beta0',Beta0,...
    'lb',[-Inf 0 -Inf]);

Method: Maximum likelihood (fmincon)
Effective Sample size:             60
Logarithmic  likelihood:     -110.477
Akaike   info criterion:      226.954
Bayesian info criterion:      233.287
           |      Coeff       Std Err    t Stat    Prob 
--------------------------------------------------------
 c(1)      |   0.59436       0.09408     6.31738  0     
 c(2)      |   1.52554       0.10758    14.17991  0     
 y <- z(1) | -24.26161       1.55730   -15.57930  0     
           |                                            
           |    Final State   Std Dev     t Stat   Prob 
 x(1)      |   2.54764        0           Inf     0     

EstMdl является моделью dssm, и можно получить доступ к ее свойствам с помощью записи через точку.

Сглаживайте предполагаемую рассеянную модель в пространстве состояний. EstMdl не хранит данные или коэффициенты регрессии, таким образом, необходимо передать в них в использовании аргументов пары "имя-значение" 'Predictors' и 'Beta', соответственно. Постройте сглаживавшие состояния.

SmoothedX = smooth(EstMdl,y,'Predictors',Z,'Beta',estParams(end));

figure
plot(dates(end-(T-1)+1:end),SmoothedX);
xlabel('Period')
ylabel('Change in the unemployment rate')
title('Smoothed Change in the Unemployment Rate')
axis tight

Оцените рассеянную модель в пространстве состояний, сглаживайте состояния, и затем извлеките другие оценки от выходного аргумента Output.

Рассмотрите рассеянную модель в пространстве состояний

Переменная состояния $$x_{1,t}$$ модель AR (1) с авторегрессивным коэффициентом $$\varphi$$. $$x_{2,t}$$ случайный обход. Воздействия $$u_{1,t}$$ и $$u_{2,t}$$ независимые Гауссовы случайные переменные со средним значением 0 и стандартными отклонениями $${\sigma }_1$$ и $${\sigma }_2$$, соответственно. Наблюдение $$y_t$$ безошибочная сумма $$x_{1,t}$$ и $$x_{2,t}$$.

Сгенерируйте данные из модели в пространстве состояний. Чтобы моделировать данные, предположите что объем выборки $$T=\mathrm{}\mathrm{}100$$, $$\varphi =0.6$$, $${\sigma }_1=0.2$$, $${\sigma }_2=0.1$$, и $$x_{1,0}=x_{2,0}=2$$.

rng(1); % For reproducibility
T = 100;
ARMdl = arima('AR',0.6,'Constant',0,'Variance',0.2^2);
x1 = simulate(ARMdl,T,'Y0',2);
u3 = 0.1*randn(T,1);
x3 = cumsum([2;u3]);
x3 = x3(2:end);
y = x1 + x3;

Задайте содействующие матрицы модели в пространстве состояний. Чтобы указать на неизвестные параметры, используйте значения NaN.

A = [NaN 0; 0 1];
B = [NaN 0; 0 NaN];
C = [1 1];

Создайте рассеянную модель в пространстве состояний, которая описывает модель выше. Задайте это $$x_{1,t}$$ и $$x_{2,t}$$ имейте рассеянные дистрибутивы начального состояния.

StateType = [2 2];
Mdl = dssm(A,B,C,'StateType',StateType);

Оцените неизвестные параметры Mdl. Выберите начальные значения параметров для оптимизации. Укажите, что стандартные отклонения ограничиваются быть положительными, но все другие параметры являются неограниченным использованием аргумента пары "имя-значение" 'lb'.

params0 = [0.01 0.1 0.01]; % Initial values chosen arbitrarily
EstMdl = estimate(Mdl,y,params0,'lb',[-Inf 0 0]);

Method: Maximum likelihood (fmincon)
Effective Sample size:             98
Logarithmic  likelihood:      3.44283
Akaike   info criterion:    -0.885655
Bayesian info criterion:      6.92986
      |     Coeff      Std Err   t Stat     Prob  
--------------------------------------------------
 c(1) | 0.54134       0.20494    2.64145  0.00826 
 c(2) | 0.18439       0.03305    5.57897   0      
 c(3) | 0.11783       0.04347    2.71039  0.00672 
      |                                           
      |  Final State   Std Dev    t Stat    Prob  
 x(1) | 0.24884       0.17168    1.44943  0.14722 
 x(2) | 1.73762       0.17168   10.12121   0      

Параметры близко к их истинным значениям.

Сглаживайте состояния EstMdl и запросите весь другой доступный вывод.

[X,logL,Output] = smooth(EstMdl,y);

X является T-by-2 матрица сглаживавших состояний, logL является оптимизированным значением логарифмической вероятности финала, и Output является массивом структур, содержащим различные оценки, что Фильтр Калмана требует. Перечислите поля output с помощью fields.

fields(Output)

ans = 9x1 cell array
    {'LogLikelihood'          }
    {'SmoothedStates'         }
    {'SmoothedStatesCov'      }
    {'SmoothedStateDisturb'   }
    {'SmoothedStateDisturbCov'}
    {'SmoothedObsInnov'       }
    {'SmoothedObsInnovCov'    }
    {'KalmanGain'             }
    {'DataUsed'               }

Преобразуйте Output в таблицу.

OutputTbl = struct2table(Output);
OutputTbl(1:10,1:4) % Display first ten rows of first four variables

ans=10×4 table
    LogLikelihood    SmoothedStates    SmoothedStatesCov    SmoothedStateDisturb
    _____________    ______________    _________________    ____________________

      []              []                 []                     []              
      []              []                 []                     []              
      [ 0.1827]       [2x1 double]       [2x2 double]           [2x1 double]    
      [ 0.0972]       [2x1 double]       [2x2 double]           [2x1 double]    
      [ 0.4472]       [2x1 double]       [2x2 double]           [2x1 double]    
      [ 0.2073]       [2x1 double]       [2x2 double]           [2x1 double]    
      [ 0.5167]       [2x1 double]       [2x2 double]           [2x1 double]    
      [ 0.2389]       [2x1 double]       [2x2 double]           [2x1 double]    
      [ 0.5064]       [2x1 double]       [2x2 double]           [2x1 double]    
      [-0.0105]       [2x1 double]       [2x2 double]           [2x1 double]    

Первые две строки таблицы содержат пустые ячейки или нули. Они соответствуют наблюдениям, требуемым инициализировать рассеянный Фильтр Калмана. Таким образом, SwitchTime равняется 2.

SwitchTime = 2;

Постройте сглаживавшие состояния и их отдельные 95% доверительных интервалов вальдового типа.

CI = nan(T,2,2);

for j = (SwitchTime + 1):T
    CovX = OutputTbl.SmoothedStatesCov{j};
    CI(j,:,1) = X(j,1) + 1.96*sqrt(CovX(1,1))*[-1 1];
    CI(j,:,2) = X(j,2) + 1.96*sqrt(CovX(2,2))*[-1 1];    
end

figure;
plot(1:T,X(:,1),'k',1:T,CI(:,:,1),'--r');
xlabel('Period');
ylabel('Smoothed states');
title('State 1 Estimates')
legend('Smoothed','95% Individual CIs');
grid on;

figure;
plot(1:T,X(:,2),'k',1:T,CI(:,:,2),'--r');
xlabel('Period');
ylabel('Smoothed states');
title('State 2 Estimates')
legend('Smoothed','95% Individual CIs');
grid on;