Предположим, что скрытый процесс является AR (1). Уравнение состояния

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

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

T = 100;
ARMdl = arima('AR',0.5,'Constant',0,'Variance',1);
x0 = 1.5;
rng(1); % For reproducibility
x = simulate(ARMdl,T,'Y0',x0);

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

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

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

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

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

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

Задайте модель в пространстве состояний с помощью содействующих матриц.

Mdl = ssm(A,B,C,D)

Mdl = 
State-space model type: ssm

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) = (0.50)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  1.33 

State types
     x1     
 Stationary 

Mdl ssm модель. Проверьте, что модель правильно задана с помощью отображения в Командном окне. Программное обеспечение выводит, что процесс состояния является стационарным. Впоследствии, программное обеспечение устанавливает среднее значение начального состояния и ковариацию к среднему значению и отклонению стационарного распределения модели AR (1).

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

где:

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

load Data_NelsonPlosser

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

isNaN = any(ismissing(DataTable),2);       % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
u = DataTable.UR(~isNaN);
T = size(gnpn,1);                          % Sample size
Z = [ones(T-1,1) diff(log(gnpn))];
y = diff(u);

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

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

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

Задайте модель в пространстве состояний с помощью ssm.

Mdl = ssm(A,B,C,D);

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

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

Method: Maximum likelihood (fmincon)
Sample size: 61
Logarithmic  likelihood:     -99.7245
Akaike   info criterion:      209.449
Bayesian info criterion:      220.003
           |      Coeff       Std Err    t Stat     Prob  
----------------------------------------------------------
 c(1)      |  -0.34098       0.29608    -1.15164  0.24948 
 c(2)      |   1.05003       0.41377     2.53771  0.01116 
 c(3)      |   0.48592       0.36790     1.32079  0.18657 
 y <- z(1) |   1.36121       0.22338     6.09358   0      
 y <- z(2) | -24.46711       1.60018   -15.29024   0      
           |                                              
           |    Final State   Std Dev     t Stat    Prob  
 x(1)      |   1.01264       0.44690     2.26592  0.02346 
 x(2)      |   0.77718       0.58917     1.31912  0.18713 

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

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

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

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

Считайте nowcasting моделью в состояниях Фильтра Независимой от времени Модели в пространстве состояний.

Сгенерируйте случайную последовательность 100 наблюдений от $$x_t$$.

T = 100;
A = 0.5;
B = 1;
C = 1;
D = 0.75;
Mdl = ssm(A,B,C,D);
rng(1); % For reproducibility
y = simulate(Mdl,T);

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

fh = 10;
yf = y((end-fh+1):end); % Holdout sample responses
y = y(1:end-fh);        % In-sample responses

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

[xhat,logL,output] = filter(Mdl,y);
xhatvar = output.FilteredStatesCov;

xhat и xhatvar 90 1 векторы из отфильтрованных состояний в выборке и соответствующих отклонений, соответственно. xhat(t) оценка $\mathit{E}\left({\mathit{x}}_{\mathit{t}}|{\mathit{y}}_1,...,{\mathit{y}}_{\mathit{t}}\right)$, и xhatvar оценка $\mathrm{Var}\left({\mathit{x}}_{\mathit{t}}|{\mathit{y}}_1,...,{\mathit{y}}_{\mathit{t}}\right)$.

Вызовите filter снова, но задайте опцию обновления в реальном времени.

[xhatRT,logLRT,outputRT] = filter(Mdl,y,'RealTimeUpdate',true);
xhatRTvar = outputRT.FilteredStatesCov;

xhatRT и xhatRTvar скаляры, представляющие оценку $\mathit{E}\left({\mathit{x}}_{90}|{\mathit{y}}_1,...,{\mathit{y}}_{90}\right)$ и его соответствующее отклонение, соответственно.

Сравните отфильтрованные состояния и отклонения периода 90.

tol = 1e-10;
areMeansEqual = (xhat(end) - xhatRT) < tol

areMeansEqual = logical
   1

areVarsEqual = (xhatvar(end) - xhatRTvar) < tol

areVarsEqual = logical
   0

areLogLsEqual = (logL - logLRT) < tol;

В последний период отфильтрованные состояния, их отклонения и логарифмическая правдоподобность равны.

Nowcast модель в горизонт прогноза путем выполнения этой процедуры в течение каждого последовательного периода:

% Initialize state and variance
state0 = xhatRT;        
var0 = xhatRTvar;

% Preallocate
xhatRTF = zeros(fh,1);
xhatRTvarF = zeros(fh,1);

for j = 1:fh
    Mdl.Mean0 = state0;
    Mdl.Cov0 = var0;
    [xhatRTF(j),~,outputRT] = filter(Mdl,yf(j),'RealTimeUpdate',true); % Alternatively, use update
    xhatRTvarF(j) = outputRT.FilteredStatesCov;
    state0 = xhatRTF(j);
    var0 = xhatRTvarF(j);
end

Отобразите на графике данные и nowcasts.

figure
plot((T-fh-20):T,[y(end-20:end); yf],'b-',(T-fh+1):T,xhatRTF,'r*-')
legend(["Data" "Nowcasts"],'Location',"best")