Оцените случайный параметр модели в пространстве состояний

В этом примере показано, как оценить случайный, авторегрессивный коэффициент состояния в модели в пространстве состояний. Таким образом, этот пример получает представление Bayesian оценки параметра модели в пространстве состояний при помощи "зигзагообразного" метода оценки.

Предположим это два состояния (x1,t и x2,t) представляйте сетевой экспорт двух стран в конце года.

  • x1,t модульный корневой процесс с отклонением воздействия σ12.

  • x2,t AR (1) процесс с неизвестным, случайным коэффициентом и отклонением воздействия σ22.

  • Наблюдение (yt) точная сумма двух сетевого экспорта. Таким образом, сетевой экспорт отдельных государств неизвестен.

Символически, истинная модель в пространстве состояний

x1,t=x1,t-1+σ1u1,tx2,t=ϕx2,t-1+σ2u2,tyt=x1,t+x2,t

Симулируйте данные

Симулируйте 100 лет сетевого экспорта от:

  • Модульный корневой процесс со средним нулем, Гауссов шумовой ряд, который имеет отклонение 0.12.

  • AR (1) процесс с авторегрессивным коэффициентом 0,6 и средний нуль, Гауссов шумовой ряд, который имеет отклонение 0.22.

  • x1,0=x2,0=0.

  • Создайте ряд наблюдения путем подведения итогов двух сетевого экспорта в год.

rng(100); % For reproducibility
T = 150;
sigma1 = 0.1;
sigma2 = 0.2;
phi = 0.6;

u1 = randn(T,1)*sigma1;
x1 = cumsum(u1);
Mdl2 = arima('AR',phi,'Variance',sigma2^2,'Constant',0);
x2 = simulate(Mdl2,T,'Y0',0);
y = x1 + x2;

figure;
plot([x1 x2 y])
legend('x_1','x_2','y','Location','Best');
ylabel('Net exports');
xlabel('Period');

Figure contains an axes. The axes contains 3 objects of type line. These objects represent x_1, x_2, y.

Зигзагообразный метод оценки

Обработка ϕ как будто это неизвестно и случайно, и используйте зигзагообразный метод, чтобы восстановить его распределение. Реализовывать зигзагообразный метод:

1. Выберите начальное значение для ϕ в интервале (-1,1), и обозначают его ϕz.

2. Создайте истинную модель в пространстве состояний, то есть, ssm объект модели, который представляет генерирующий данные процесс.

3. Используйте более сглаженную симуляцию (simsmooth) чертить случайный путь от распределения вторых сглаживавших состояний. Символически, x2,z,tP(x2,t|yt,x1,t,ϕ=ϕz).

4. Создайте другую модель в пространстве состояний, которая имеет эту форму

ϕz,t=ϕz,t-1x2,z,t=x2,z,t-1ϕz,t+0.8u2,t

В словах, ϕz,t статическое состояние и x2,z,t "наблюдаемый" ряд со временем различный коэффициент Ct=x2,z,t-1.

5. Используйте симуляцию, более сглаженную, чтобы чертить случайный путь от распределения сглаживавшего ϕz,t ряд. Символически, ϕz,tP(ϕ|x2,z,t), где x2,z,t охватывает структуру истинной модели в пространстве состояний и наблюдений. ϕz,t является статическим, таким образом, можно только зарезервировать одно значение (ϕz).

6. Повторяйте шаги 2 - 5 много раз и хранилище ϕz каждая итерация.

7. Выполните диагностические проверки на ряде симуляции. Таким образом, создайте:

  • Проследите графики определить запись в период и смешивается ли Цепь Маркова хорошо.

  • Автокорреляция строит, чтобы определить, сколько чертит удаление потребности, чтобы получить хорошо смешанную Цепь Маркова.

8. Остающийся ряд представляет, чертит от апостериорного распределения ϕ. Можно вычислить описательную статистику или построить гистограмму, чтобы определить качества распределения.

Оцените случайный коэффициент Используя зигзагообразный метод

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

phi0 = -0.3;             % Initial value of phi
Z = 1000;                % Number of times to iterate the zig-zag method
phiz = [phi0; nan(Z,1)]; % Preallocate

A = [1 0; 0 NaN];
B = [sigma1; sigma2];
C = [1 1];
Mdl = ssm(A,B,C,'StateType',[2; 0]);

Mdl является ssm объект модели. NaN действия как заполнитель для ϕ.

Выполните итерации шагов 2 - 5 зигзагообразного метода.

for j = 2:(Z + 1);
    % Draw a random path from smoothed x_2 series.
    xz = simsmooth(Mdl,y,'Params',phiz(j-1));
    % The second column of xz is a draw from the posterior distribution of x_2.
    
    % Create the intermediate state-space model.
    Az = 1;
    Bz = 0;
    Cz = num2cell(xz((1:(end - 1)),2));
    Dz = sigma2;
    Mdlz = ssm(Az,Bz,Cz,Dz,'StateType',2);
    
    % Draw a random path from the smoothed phiz series.
    phizvec = simsmooth(Mdlz,xz(2:end,2));
    phiz(j) = phizvec(1);
    % phiz(j) is a draw from the posterior distribution of phi
end

phiz Цепь Маркова. Прежде, чем анализировать апостериорное распределение ϕ, необходимо оценить, наложить ли электротермотренировку или серьезность автокорреляции в цепи.

Определите качество симуляции

Постройте график трассировки для первых 100, 500, и все случайные ничьи.

vec = [100 500 Z];
figure;
for j = 1:3;
    subplot(3,1,j)
    plot(phiz(1:vec(j)));
    title('Trace Plot for \phi');
    xlabel('Simulation number');
    axis tight;
end

Figure contains 3 axes. Axes 1 with title Trace Plot for \phi contains an object of type line. Axes 2 with title Trace Plot for \phi contains an object of type line. Axes 3 with title Trace Plot for \phi contains an object of type line.

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

График автокорреляционная функция ряда после удаления первых 20 чертит.

burnOut = 21:Z;

figure;
autocorr(phiz(burnOut));

Figure contains an axes. The axes with title Sample Autocorrelation Function contains 4 objects of type stem, line.

Автокорреляционная функция вымирает скорее быстро. Не кажется, что автокорреляция в цепи является проблемой.

Определите качества апостериорного распределения ϕ вычислительной статистикой симуляции и путем графического вывода гистограммы уменьшаемого набора случайных ничьих.

xbar = mean(phiz(burnOut))
xbar = 0.5104
xstd = std(phiz(burnOut))
xstd = 0.0988
ci = norminv([0.025,0.975],xbar,xstd); % 95% confidence interval

figure;
histogram(phiz(burnOut),'Normalization','pdf');
h = gca;
hold on;
simX = linspace(h.XLim(1),h.XLim(2),100);
simPDF = normpdf(simX,xbar,xstd); 
plot(simX,simPDF,'k','LineWidth',2);
h1 = plot([xbar xbar],h.YLim,'r','LineWidth',2);
h2 = plot([0.6 0.6],h.YLim,'g','LineWidth',2);
h3 = plot([ci(1) ci(1)],h.YLim,'--r',...
    [ci(2) ci(2)],h.YLim,'--r','LineWidth',2);
legend([h1 h2 h3(1)],{'Simulation Mean','True Mean','95% CI'});
h.XTick = sort([h.XTick xbar]);
h.XTickLabel{h.XTick == xbar} = xbar;
h.XTickLabelRotation = 90;

Figure contains an axes. The axes contains 6 objects of type histogram, line. These objects represent Simulation Mean, True Mean, 95% CI.

Апостериорное распределение ϕ примерно нормально со средним и стандартным отклонением приблизительно 0,51 и 0.1, соответственно. Истинное среднее значение ϕ 0.6, и это меньше одного стандартного отклонения справа от среднего значения симуляции.

Вычислите оценку наибольшего правдоподобия ϕ. Таким образом, обработка ϕ как фиксированный, но неизвестный параметр, и затем оценивают Mdl использование Фильтра Калмана и наибольшего правдоподобия.

[~,estParams] = estimate(Mdl,y,phi0)
Method: Maximum likelihood (fminunc)
Sample size: 150
Logarithmic  likelihood:     -10.1434
Akaike   info criterion:      22.2868
Bayesian info criterion:      25.2974
      |     Coeff       Std Err     t Stat       Prob  
-------------------------------------------------------
 c(1) |  0.53590       0.19183    2.79360      0.00521 
      |                                                
      |   Final State   Std Dev    t Stat        Prob  
 x(1) | -0.85059       0.00000   -6.45811e+08   0      
 x(2) |  0.00454        0          Inf          0      
estParams = 0.5359

MLE ϕ 0.54. Обе оценки в одном стандартном отклонении или стандартной погрешности от истинного значения ϕ.

Смотрите также

| | | |

Связанные примеры

Больше о

Для просмотра документации необходимо авторизоваться на сайте