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

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

Предположим это два состояния ( $x_{1, t}$ и $x_{2, t}$ ) представляйте сетевой экспорт двух стран в конце года.

$x_{1, t}$ модульный корневой процесс с отклонением воздействия $σ_{1}^{2}$ .
$x_{2, t}$ AR (1) процесс с неизвестным, случайным коэффициентом и отклонением воздействия $σ_{2}^{2}$ .
Наблюдение ( $y_{t}$ ) точная сумма двух сетевого экспорта. Таким образом, сетевой экспорт отдельных государств неизвестен.

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

$\begin{array}{l} x_{1, t} = x_{1, t - 1} + σ_{1} u_{1, t} \\ x_{2, t} = ϕ x_{2, t - 1} + σ_{2} u_{2, t} \\ y_{t} = x_{1, t} + x_{2, t} \end{array}$

Моделируйте данные

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

Модульный корневой процесс со средним нулем, Гауссов шумовой ряд, который имеет отклонение $0 . 1^{2}$ .
AR (1) процесс с авторегрессивным коэффициентом 0,6 и средний нуль, Гауссов шумовой ряд, который имеет отклонение $0 . 2^{2}$ .
$x_{1, 0} = x_{2, 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');

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

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

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

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

3. Используйте симуляцию, более сглаженную (simsmooth), чтобы чертить случайный путь от распределения вторых сглаживавших состояний. Символически, $x_{2, z, t} \sim P (x_{2, t} | y_{t}, x_{1, t}, ϕ = ϕ_{z})$ .

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

$\begin{array}{l} ϕ_{z, t} = ϕ_{z, t - 1} \\ x_{2, z, t} = x_{2, z, t - 1} ϕ_{z, t} + 0.8 u_{2, t} \end{array}$

В словах, $ϕ_{z, t}$ статическое состояние и $x_{2, z, t}$ "наблюдаемый" ряд со временем переменный коэффициент $C_{t} = x_{2, z, t - 1}$ .

5. Используйте симуляцию, более сглаженную, чтобы чертить случайный путь от распределения сглаживавшего $ϕ_{z, t}$ ряд. Символически, $ϕ_{z, t} \sim P (ϕ | x_{2, z, t})$ , где $x_{2, 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

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

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

burnOut = 21:Z;

figure;
autocorr(phiz(burnOut));

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

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

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;

Апостериорное распределение $ϕ$ примерно нормально со средним и стандартным отклонением приблизительно 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. Обе оценки в одном стандартном отклонении или стандартной погрешности от истинного значения $ϕ$ .

Документация