Matlab马尔可夫区制转换动态回归模型估计GDP增长率

最近我们被要求撰写关于马尔可夫区制转换动态回归模型的研究报告，包括一些图形和统计输出。

相关视频：马尔可夫链原理可视化解释与R语言区制转换Markov regime switching实例

马尔可夫链原理可视化解释与R语言区制转换Markov regime switching实例

，时长07:25

相关视频

马尔可夫链蒙特卡罗方法MCMC原理与R语言实现

，时长08:47

本文估计实际GDP增长率的两状态Markov区制转换动态回归模型 。

创建模型进行估计

通过指定转移矩阵和两个区制的AR（0）（仅常数）子模型的两状态离散时间马尔可夫链，为朴素估计量创建马尔可夫转换动态回归模型。标记状态。

mc(NaN(2),'StateNames',["增长" "衰退"]);

加载和预处理数据

加载GDP数据集。

Data 包含1947：Q1-2005：Q2期间实际GDP的季度数据。估计周期  为1947：Q2-2004：Q2。

通过以下方式将数据转换为年度序列：

1. 在估计期内将数据转换为季度比率

2. 将季度比率年度化

diff(Data(2:230))./Data(2:229); % 季度比率
100*((1 + qrate).^4 - 1);       % 年度比率

估计模型

模型拟合Mdl 年利率序列 arate。

 estimate(Mdl,Mdl0,arate);

EstMdl 是估计的（完全指定的）马尔可夫转换动态回归模型。 EstMdl.Switch 是估计的离散时间马尔可夫链模型（dtmc 对象）， EstMdl.Submodels 是估计的单变量VAR（0）模型（varm 对象）的向量。

显示估计的特定于状态的动态模型。

  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

  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

显示估计的状态转移矩阵。

2×2

    0.9088    0.0912
    0.2303    0.7697

估计后的EM算法

在估计马尔可夫转换动态回归模型中考虑模型和数据 。

创建部分指定模型进行估计。

创建完全指定的模型，其中包含估计过程的初始参数值。

加载并预处理数据。

diff(Data(2:230))./Data(2:229); 
100*((1 + qrate).^4 - 1);       

使模型拟合数据。当估计过程终止时，绘制对数似然比与迭代步骤。

Plot(Mdl,Mdl0);

使模型拟合模拟数据

使用来自已知数据生成过程（DGP）的模拟数据评估估计准确性。本示例使用任意参数值。

为DGP创建模型

为转换区制创建一个完全指定的两状态离散时间马尔可夫链模型。

P = [0.7 0.3; 0.1 0.9];

对于每个状态，为过程创建一个完全指定的AR（1）模型。

% 常数
C1 = 4;
C2 = -1;

% 自回归系数
AR1 = 0.5;
AR2 = 0.3;

% 方差
V1 = 3;
V2 = 2;

% AR 子模型
arima('Constant',C1,'AR',AR1,'Variance',V1);

为DGP创建完全指定的Markov转换动态回归模型。

模拟来自DGP的路径

从DGP生成10条长度为1000的随机路径。

rng(1); % 重现性
N = 10;
n = 1000;
simulate(DGP,n,'Numpaths',N);

Data 是模拟的1000 x 10矩阵。

创建估计模型

创建一个部分指定的马尔可夫转换动态回归模型，该模型具有与数据生成过程相同的结构，但是指定了未知的转移矩阵和未知的子模型系数。

创建包含初始值的模型

创建一个完全指定的马尔可夫转换动态回归模型，该模型具有与相同的结构 Mdl，但是将所有可估计的参数都设置为初始值。

P0 = 0.5*ones(2);
dtmc(P0);

ms(mc0,[mdl01,mdl02]);

估计模型

使模型拟合每个模拟路径。对于每条路径，在EM算法的每次迭代中绘制对数似然图。

figure
hold on

for i = 1:N

    estimate(Data(:,i),'Plot',true);
    
    
end

hold off

评估准确性

计算每个估计参数的蒙特卡洛平均值。

将总体参数与相应的蒙特卡洛估计进行比较。

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 = 2×2

    0.7000    0.3000
    0.1000    0.9000

PEstimate = 2×2

    0.7065    0.2935
    0.1023    0.8977

预采样数据

考虑 估计马尔可夫转换动态回归模型中的数据，但假设关注期间为1960：Q1–2004：Q2。另外，考虑向每个子模型添加一个自回归项。

创建部分指定的马尔可夫转换动态回归模型进行估计。指定AR（1）子模型。

arima(1,0,0);
ms(mc,[mdl; mdl]);

由于子模型是AR（1），因此每个子模型都需要进行一次预采样观察以初始化其动态分量以进行估计。

创建包含用于估计过程的初始参数值的模型。

P0 = 0.5*ones(2);
mc(P0,'StateNames');

加载数据。将整个集合转换为年化利率序列。

使用与年率序列相关的日期来确定预采样和估计采样周期。由于转换应用了一阶差分，因此必须从原始样本中删除第一个观察日期。

dates = datetime(dates(2:end),'ConvertFrom','datenum',...
    'Format','yyyy:QQQ','Locale','en_US');
estPrd = datetime(["1960:Q2" "2004:Q2"],'InputFormat','yyyy:QQQ',...
    'Format','yyyy:QQQ','Locale','en_US');

使模型拟合估计样本数据。指定预采样观测值，并在估计过程终止时在每次迭代时绘制对数似然度。

访问预期的平滑状态概率和对数似然

在估计马尔可夫转换动态回归模型中考虑模型和数据 。

创建部分指定的模型进行估计。

创建完全指定的模型，其中包含估计过程的初始参数值。

加载并预处理数据。

使模型拟合数据。当算法终止时，返回预期的平滑状态概率和对数似然。

[EstMdl,SS,logL] = estimate(Mdl,Mdl0,arate);

SS 是预期平滑状态概率的228 x 2矩阵；行对应于估计样本中的周期，列对应于方案。 logL 是最终的对数似然。

显示估计样本中最后一个时期的预期平滑状态概率，并显示最终对数似然。

ans = 1×2

    0.8985    0.1015

logL = -639.4962

执行约束估计

将模拟数据拟合到具有VARX子模型的Markov转换动态回归模型。指定用于估计的相等约束。

为DGP创建模型

为转换区制创建一个完全指定的三态离散时间马尔可夫链模型。

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

对于每种状态，为响应过程创建一个完全指定的VARX（1）模型。为所有子模型指定相同的模型常数和滞后1 AR系数矩阵。对于每个模型，为一个外生变量指定不同的回归系数。

% 常数
C = [1;-1];

% 自回归系数
AR = {[0.6 0.1; 0.4 0.2]};

% 回归系数
Beta1 = [0.2;-0.4];


% VAR 子模型
dgp = varm('Constant',C,'AR',AR,'Covariance',5*eye(2));

为DGP创建完全指定的Markov转换动态回归模型。

ms(mc,[dgp1; dgp2; dgp3]);

模拟来自DGP的数据

通过从均值0和方差100的高斯分布中生成1000个观测值来模拟外生序列的数据。

rng(1); % 重现性
X = 10*randn(1000,1);

从DGP生成长度为1000的随机路径。为子模型回归指定模拟的外部数据。

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

Data 是模拟的1000 x 1向量。

创建估计模型

创建一个部分指定的马尔可夫转换动态回归模型，该模型具有与数据生成过程相同的结构，但是指定了未知的转换矩阵和未知的回归系数。指定常数和AR系数矩阵的真值。

ms(mcEst,[mdl; mdl; mdl]);

由于常数和AR系数矩阵的值被指定在 Mdl， estimate 将它们作为用于估计等式约束。

创建包含初始值的模型

创建具有与相同结构的完全指定的马尔可夫转换动态回归模型 Mdl，但将所有可估计参数设置为初始值，并将具有相等约束的参数设置为中指定的值 Mdl。

估计模型

使模型拟合模拟数据。指定回归的外部数据。在EM算法的每次迭代中绘制对数似然。

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

评估准确性

将估计的回归系数向量和转移矩阵与其真实值进行比较。

Beta1 = 2×1

    0.2000
   -0.4000

Beta1Estimate = 2×1

    0.1596
   -0.4040

Beta2 = 2×1

    0.6000
   -1.0000

Beta2Estimate = 2×1

    0.5888
   -0.9771

Beta3 = 2×1

    0.9000
   -1.3000

Beta3Estimate = 2×1

    0.8987
   -1.2991

P = 3×3

    0.8000    0.1000    0.1000
    0.2000    0.6000    0.2000
         0    0.1000    0.9000

PEstimate = 3×3

    0.7787    0.0856    0.1357
    0.1366    0.6906    0.1727
    0.0086    0.0787    0.9127

---

最受欢迎的见解

1.用R语言模拟混合制排队随机服务排队系统

2.R语言中使用排队论预测等待时间

3.R语言中实现马尔可夫链蒙特卡罗MCMC模型

4.R语言中的马尔科夫机制转换(Markov regime switching)模型

5.matlab贝叶斯隐马尔可夫hmm模型

6.用R语言模拟混合制排队随机服务排队系统

7.Python基于粒子群优化的投资组合优化

8.R语言马尔可夫转换模型研究交通伤亡人数事故预测

9.用机器学习识别不断变化的股市状况——隐马尔可夫模型的应用