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

Симулируйте 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$$ как будто это неизвестно и случайно, и используйте зигзагообразный метод, чтобы восстановить его распределение. Реализовывать зигзагообразный метод:

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