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

где $${}_{\mathrm{ut}}$$ - гауссов со средним значением 0 и стандартным отклонением 1.

Создайте случайную серию из 100 наблюдений из $${}_{\mathrm{xt}}$$, предполагая, что серия начинается с 1.5.

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

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

где $${}_{\mathrm{\alpha 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). Символично и в форме state-space модель

где:

Загрузите набор данных Нельсона-Плоссера, который содержит, среди прочего, уровень безработицы и ряды 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' аргументы пары имя-значение соответственно. Просмотрите оценки и всю диагностическую информацию по оптимизации. Ограничьте оценку $$\lambda$$ всеми положительными, вещественными числами.

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 выходной аргумент.

Рассмотрим модель диффузного состояния и пространства

Переменная состояния $${}_{\mathrm{x1},}$$t представляет собой модель AR (1) с авторегрессионным $$\mathrm{коэффициентом}$$$${}_{\mathrm{}\mathrm{x2}}$$, t - случайная прогулка. $${}_{\mathrm{Возмущения}}$$u1, t $${}_{\mathrm{}и}$$u2, t - независимые гауссовы случайные величины со средними 0 и стандартными $${}_{\mathrm{}}$$отклонениями start1 $${}_{\mathrm{}}$$и start2 соответственно. $${}_{}$$Наблюдение yt представляет собой безошибочную $${\mathrm{сумму}}_{\mathrm{}}$$x1, $$t_{\mathrm{}}$$и x2, t.

Создание данных из модели пространства состояний. Чтобы смоделировать данные, предположим, что размер выборки $$T\mathrm{=}\mathrm{}\mathrm{}$$100$$\mathrm{}\mathrm{}$$$${}_{\mathrm{}}\mathrm{}\mathrm{}$$$${}_{\mathrm{}}\mathrm{}\mathrm{}$$$${}_{\mathrm{}\mathrm{}}{}_{\mathrm{}\mathrm{}}\mathrm{}$$,

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;

Задайте матрицы коэффициентов модели state-space. Для указания неизвестных параметров используйте NaN значения.

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

Создайте модель диффузного состояния-пространства, описывающую модель выше. Укажите, что $${}_{\mathrm{x1},}$$t и $${}_{\mathrm{}x2,}$$ 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;