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

Смоделировать 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');

Метод оценки Zig-Zag

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

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

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

3. Используйте сглаживание моделирования (simsmooth), чтобы нарисовать случайный путь из распределения вторых сглаженных состояний. Символически, $${}_{\mathrm{}}\mathrm{x2,z,t\sim P}{(}_{\mathrm{}x2,}{}_{}{\mathrm{t\; \text{'}yt}}_{\mathrm{,}}x1,t_{}$$ ,

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

В словах, $${}_{\mathrm{\delta \; z},}$$t - статическое состояние, а $${}_{\mathrm{}x2,z}$$, t - «наблюдаемый» ряд с изменяющимся во времени $${}_{\mathrm{коэффициентом}}{}_{\mathrm{Ct}=x2\mathrm{,}}$$z, t-1.

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

6. Повторите шаги 2-5 много раз и сохраните (сохраняйте, сохраняйте) $${}_{\mathrm{для}}$$ каждой итерации.

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

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

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

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

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 выступает в качестве местозаполнителя $$\mathrm{для}$$

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

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;

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

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