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

Создайте модель динамической регрессии с переключением Маркова для наивного оценщика путем определения дискретной цепи Маркова двух состояний с неизвестной матрицей перехода и подмоделей 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]);

Потому что значения констант и матриц коэффициентов заданы в 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