Создание частично заданной модели для оценки

Создайте модель динамической регрессии с коммутацией Маркова для наивного оценщика, задав дискретную цепочку Маркова с двумя состояниями с неизвестной матрицей перехода и подмоделями AR (0) (только константы) для обоих режимов. Маркировать режимы.

P = NaN(2);
mc = dtmc(P,'StateNames',["Expansion" "Recession"]);
mdl = arima(0,0,0);
Mdl = msVAR(mc,[mdl; mdl]);

Mdl является частично указанным msVAR объект. NaN-значимые элементы Switch и SubModels свойства указывают оцениваемые параметры.

Создание полностью заданной модели, содержащей начальные значения

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

P0 = 0.5*ones(2);
mc0 = dtmc(P0,'StateNames',Mdl.StateNames);
mdl01 = arima('Constant',1,'Variance',1);
mdl02 = arima('Constant',-1,'Variance',1);
Mdl0 = msVAR(mc0,[mdl01; mdl02]);

Mdl0 является полностью указанным msVAR объект.

Загрузка и предварительная обработка данных

Загрузите набор данных ВВП США.

load Data_GDP

Data содержит квартальные измерения реального ВВП США в период 1947: Q1-2005: Q2. Период оценки в [1] равен 1947:Q2-2004:Q2. Для получения дополнительной информации о наборе данных введитеDescription в командной строке.

Преобразуйте данные в годовой ряд скоростей, выполнив следующие действия:

qrate = diff(Data(2:230))./Data(2:229); % Quarterly rate
arate = 100*((1 + qrate).^4 - 1);       % Annualized rate

Модель оценки

Подгонка модели Mdl в ряды годовых ставок arate. Определить Mdl0 в качестве модели, содержащей начальные оценочные значения параметров.

EstMdl = estimate(Mdl,Mdl0,arate);

EstMdl - оценочная (полностью заданная) модель динамической регрессии с коммутацией Маркова. EstMdl.Switch - оценочная дискретно-временная модель цепи Маркова (dtmc объект), и EstMdl.Submodels является вектором оцененных одномерных моделей VAR (0) (varm объекты).

Отображение расчетных динамических моделей, зависящих от состояния.

EstMdlExp = EstMdl.Submodels(1)

EstMdlExp = 
  varm with properties:

     Description: "ARIMA(0,0,0) Model (Gaussian Distribution)"
     SeriesNames: "Y1" 
       NumSeries: 1
               P: 0
        Constant: 4.90146
              AR: {}
           Trend: 0
            Beta: [1×0 matrix]
      Covariance: 12.087
 
   1-Dimensional VAR(0) Model

EstMdlRec = EstMdl.Submodels(2)

EstMdlRec = 
  varm with properties:

     Description: "ARIMA(0,0,0) Model (Gaussian Distribution)"
     SeriesNames: "Y1" 
       NumSeries: 1
               P: 0
        Constant: 0.0084884
              AR: {}
           Trend: 0
            Beta: [1×0 matrix]
      Covariance: 12.6876
 
   1-Dimensional VAR(0) Model

Отображение оценочной матрицы перехода состояния.

EstP = EstMdl.Switch.P

EstP = 2×2

    0.9088    0.0912
    0.2303    0.7697

Создание модели для DGP

Создайте полностью заданную дискретную модель цепи Маркова с двумя состояниями для механизма переключения.

P = [0.7 0.3; 0.1 0.9];
mc = dtmc(P);

Для каждого состояния создайте полностью заданную модель AR (1) для процесса ответа.

% Constants
C1 = 5;
C2 = -2;

% Autoregression coefficients
AR1 = 0.4;
AR2 = 0.2;

% Variances
V1 = 4;
V2 = 2;

% AR Submodels
dgp1 = arima('Constant',C1,'AR',AR1,'Variance',V1);
dgp2 = arima('Constant',C2,'AR',AR2,'Variance',V2);

Создайте полностью заданную модель динамической регрессии с коммутацией Маркова для DGP.

DGP = msVAR(mc,[dgp1,dgp2]);

Моделирование путей ответа от DGP

Сгенерируйте 10 случайных путей отклика длиной 1000 из DGP.

rng(1); % For reproducibility
N = 10;
n = 1000;
Data = simulate(DGP,n,'Numpaths',N);

Data представляет собой матрицу смоделированных откликов 1000 на 10.

Создание модели для оценки

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

PEst = NaN(2);
mcEst = dtmc(PEst);
mdl = arima(1,0,0);
Mdl = msVAR(mcEst,[mdl; mdl]);

Создание модели, содержащей начальные значения

Создание полностью заданной модели динамической регрессии с коммутацией Маркова, которая имеет ту же структуру, что и Mdl, но установите для всех оцениваемых параметров начальные значения.

P0 = 0.5*ones(2);
mc0 = dtmc(P0);
mdl01 = arima('Constant',1,'AR',0.5,'Variance',2);
mdl02 = arima('Constant',-1,'AR',0.5,'Variance',1);
Mdl0 = msVAR(mc0,[mdl01,mdl02]);

Модели оценки

Подберите модель для каждого моделируемого пути. Для каждого пути постройте график логической привязки на каждой итерации алгоритма EM.

c1 = zeros(N,1);
c2 = zeros(N,1);
v1 = zeros(N,1);
v2 = zeros(N,1);
ar1 = zeros(N,1);
ar2 = zeros(N,1);
PStack = zeros(2,2,N);

figure
hold on

for i = 1:N

    EstModel = estimate(Mdl,Mdl0,Data(:,i),'IterationPlot',true);
    
    c1(i) = EstModel.Submodels(1).Constant;
    c2(i) = EstModel.Submodels(2).Constant;
    v1(i) = EstModel.Submodels(1).Covariance;
    v2(i) = EstModel.Submodels(2).Covariance;
    ar1(i) = EstModel.Submodels(1).AR{1};
    ar2(i) = EstModel.Submodels(2).AR{1};
    PStack(:,:,i) = EstModel.Switch.P;
    
end

hold off

Оценка точности

Вычислите среднее значение Монте-Карло для каждого расчетного параметра.

c1Mean = mean(c1);
c2Mean = mean(c2);
v1Mean = mean(v1);
v2Mean = mean(v2);
ar1Mean = mean(ar1);
ar2Mean = mean(ar2);
PMean = mean(PStack,3);

Сравните параметры населения с соответствующими оценками Монте-Карло.

DGPvsEstimate = [...
    C1  c1Mean
    C2  c2Mean
    V1  v1Mean
    V2  v2Mean
    AR1 ar1Mean
    AR2 ar2Mean]

DGPvsEstimate = 6×2

    5.0000    5.0260
   -2.0000   -1.9615
    4.0000    3.9710
    2.0000    1.9903
    0.4000    0.4061
    0.2000    0.2017

P

P = 2×2

    0.7000    0.3000
    0.1000    0.9000

PEstimate = PMean

PEstimate = 2×2

    0.7065    0.2935
    0.1023    0.8977

Создание модели для DGP

Создайте полностью заданную модель цепи Маркова с тремя состояниями дискретного времени для механизма переключения.

P = [0.8 0.1 0.1; 0.2 0.6 0.2; 0 0.1 0.9];
mc = dtmc(P);

Для каждого состояния создайте полностью заданную модель VARX (1) для процесса ответа. Укажите одну и ту же константу модели и матрицу коэффициентов AR с запаздыванием 1 для всех подмоделей. Для каждой модели укажите отдельный коэффициент регрессии 2 на 1 для одной экзогенной переменной.

% Constants
C = [1;-1];

% Autoregression coefficients
AR = {[0.6 0.1; 0.4 0.2]};

% Regression coefficients
Beta1 = [0.2;-0.4];
Beta2 = [0.6;-1.0];
Beta3 = [0.9;-1.3];

% VAR Submodels
dgp = varm('Constant',C,'AR',AR,'Covariance',5*eye(2));
dgp1 = dgp;
dgp1.Beta = Beta1;
dgp2 = dgp;
dgp2.Beta = Beta2;
dgp3 = dgp;
dgp3.Beta = Beta3;

Создайте полностью заданную модель динамической регрессии с коммутацией Маркова для DGP.

DGP = msVAR(mc,[dgp1; dgp2; dgp3]);

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

Моделирование данных для экзогенного ряда путем генерации 1000 наблюдений из гауссова распределения со средним значением 0 и дисперсией 100.

rng(1); % For reproducibility
X = 10*randn(1000,1);

Генерировать случайный путь отклика длиной 1000 из DGP. Укажите смоделированные экзогенные данные для компонентов регрессии подмодели.

Data = simulate(DGP,1000,'X',X);

Data является вектором смоделированных откликов 1000 на 1.

Создание модели для оценки

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

PEst = NaN(3);
mcEst = dtmc(PEst);
mdl = varm(2,1);
mdl.Constant = C;
mdl.AR = AR;
mdl.Beta = NaN(2,1);
Mdl = msVAR(mcEst,[mdl; mdl; mdl]);

Поскольку значения констант и матриц коэффициентов AR указаны в Mdl, estimate рассматривает их как ограничения равенства для оценки.

Создание модели, содержащей начальные значения

Создание полностью заданной модели динамической регрессии с коммутацией Маркова, которая имеет ту же структуру, что и Mdl, но установить для всех оцениваемых параметров начальные значения и установить параметры с ограничениями равенства их значениям, указанным в Mdl.

P0 = (1/3)*ones(3);
mc0 = dtmc(P0);
mdl0 = varm('Constant',C,'AR',AR,'Covariance',eye(2));
mdl01 = mdl0;
mdl01.Beta = [0.1;-0.1];
mdl02 = mdl0;
mdl02.Beta = [0.5;-0.5];
mdl03 = mdl0;
mdl03.Beta = [1;-1];
Mdl0 = msVAR(mc0,[mdl01; mdl02; mdl03]);

Модель оценки

Подгонка модели к моделируемым данным. Укажите экзогенные данные для компонента регрессии. Постройте график учета на каждой итерации алгоритма EM.

figure
EstMdl = estimate(Mdl,Mdl0,Data,'X',X,'IterationPlot',true);

Оценка точности

Сравните оцененные векторы коэффициентов регрессии и матрицу перехода с их истинными значениями.

Beta1

Beta1 = 2×1

    0.2000
   -0.4000

Beta1Estimate = EstMdl.Submodels(1).Beta

Beta1Estimate = 2×1

    0.1596
   -0.4040

Beta2

Beta2 = 2×1

    0.6000
   -1.0000

Beta2Estimate = EstMdl.Submodels(2).Beta

Beta2Estimate = 2×1

    0.5888
   -0.9771

Beta3

Beta3 = 2×1

    0.9000
   -1.3000

Beta3Estimate = EstMdl.Submodels(3).Beta

Beta3Estimate = 2×1

    0.8987
   -1.2991

P

P = 3×3

    0.8000    0.1000    0.1000
    0.2000    0.6000    0.2000
         0    0.1000    0.9000

PEstimate = EstMdl.Switch.P

PEstimate = 3×3

    0.7787    0.0856    0.1357
    0.1366    0.6906    0.1727
    0.0086    0.0787    0.9127