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

где $$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. Постройте истинные значения состояния и отфильтрованные оценки состояния.

filteredX = filter(Mdl,y);

figure
plot(1:T,x,'-k',1:T,filteredX,':r','LineWidth',2)
title({'State Values'})
xlabel('Period')
ylabel('State')
legend({'True state values','Filtered 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 ssm модель, и можно получить доступ к ее свойствам с помощью записи через точку.

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

filteredX = filter(EstMdl,y,'Predictors',Z,'Beta',estParams(end));

figure
plot(dates(end-(T-1)+1:end),filteredX);
xlabel('Period')
ylabel('Change in the unemployment rate')
title('Filtered 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=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] = filter(EstMdl,y);

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

fields(Output)

ans = 9x1 cell
    {'LogLikelihood'      }
    {'FilteredStates'     }
    {'FilteredStatesCov'  }
    {'ForecastedStates'   }
    {'ForecastedStatesCov'}
    {'ForecastedObs'      }
    {'ForecastedObsCov'   }
    {'KalmanGain'         }
    {'DataUsed'           }

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

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

ans=10×5 table
    LogLikelihood    FilteredStates    FilteredStatesCov    ForecastedStates    ForecastedStatesCov
    _____________    ______________    _________________    ________________    ___________________

    {0x0 double}      {0x0 double}       {0x0 double}         {0x0 double}         {0x0 double}    
    {0x0 double}      {0x0 double}       {0x0 double}         {0x0 double}         {0x0 double}    
    {[  0.1827]}      {2x1 double}       {2x2 double}         {2x1 double}         {2x2 double}    
    {[  0.0972]}      {2x1 double}       {2x2 double}         {2x1 double}         {2x2 double}    
    {[  0.4472]}      {2x1 double}       {2x2 double}         {2x1 double}         {2x2 double}    
    {[  0.2073]}      {2x1 double}       {2x2 double}         {2x1 double}         {2x2 double}    
    {[  0.5167]}      {2x1 double}       {2x2 double}         {2x1 double}         {2x2 double}    
    {[  0.2389]}      {2x1 double}       {2x2 double}         {2x1 double}         {2x2 double}    
    {[  0.5064]}      {2x1 double}       {2x2 double}         {2x1 double}         {2x2 double}    
    {[ -0.0105]}      {2x1 double}       {2x2 double}         {2x1 double}         {2x2 double}    

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

SwitchTime = 2;

Постройте отфильтрованные и предсказанные состояния.

ForeX = cell2mat(OutputTbl.ForecastedStates')'; % Orient forecasted states
ForeX = [zeros(SwitchTime,2);ForeX]; % Include zeros for initialization

figure;
plot(1:T,X(:,1),'r',1:T,ForeX(:,1),'b');
xlabel('Period');
ylabel('State estimate');
title('State 1 Estimates')
legend('Filtered','Forecasted');
grid on;

figure;
plot(1:T,X(:,2),'r',1:T,ForeX(:,2),'b');
xlabel('Period');
ylabel('State estimate');
title('State 2 Estimates')
legend('Filtered','Forecasted');
grid on;