Моделирование данных

Моделирование 100-летних чистых экспортов из:

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');

Метод оценки Зиг-Зага

Обработка $$\varphi$$ как будто он неизвестен и случайен, и используйте метод zig-zag, чтобы восстановить его распределение. Для реализации метода zig-zag:

1. Выберите начальное значение для $$\varphi$$ в интервале (-1,1), и обозначить его $${\varphi }_z$$.

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

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

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

На словах, $${\varphi }_{z,t}$$ является статическим состоянием и $$x_{2,z,t}$$ - «наблюдаемый» ряд с изменяющимся во времени коэффициентом $$C_t=x_{2,z,t-1}$$.

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

6. Повторите шаги 2-5 много раз и храните $${\varphi }_z$$ каждую итерацию.

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

8. Оставшаяся серия представляет черты из апостериорного распределения $$\varphi$$. Можно вычислить описательную статистику или построить гистограмму, чтобы определить качества распределения.

Оценка случайного коэффициента с использованием метода Зиг-Зага

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

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 объект модели. The NaN действует как заполнитель для $$\varphi$$.

Повторите шаги 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 является марковской цепью. Перед анализом апостериорного распределения $$\varphi$$необходимо оценить, вводить ли период ожога или тяжесть автокорреляции в цепи.

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

Нарисуйте график трассировки для первых 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));

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

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

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;

Апостериорное распределение $$\varphi$$ является примерно нормальным со средним и стандартным отклонениями приблизительно 0,51 и 0,1, соответственно. Истинное среднее значение $$\varphi$$ равен 0,6, и это меньше одного стандартного отклонения справа от среднего значения симуляции.

Вычислите максимальную оценку правдоподобия $$\varphi$$. То есть лечить $$\varphi$$ как фиксированный, но неизвестный параметр, и затем оцените 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 of $$\varphi$$ составляет 0,54. Обе оценки находятся в пределах одного стандартного отклонения или стандартной ошибки от истинного значения $$\varphi$$.