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

где $$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'})

Истинные значения и оценки фильтра примерно одинаковы, за исключением первого отфильтрованного состояния, которое равняется нулю.

Предположим, что интерес представляет линейная связь между уровнем безработицы и номинальным валовым национальным продуктом (ННП). Предположим далее, что уровень безработицы является серией 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-by-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;