TVP-VAR模型的实证分析与应用

一、TVP-VAR模型简介

TVP-VAR（Time Varying Parameter-Stochastic Volatility-Vector Auto Regression）是一种扩展的向量自回归模型，与传统VAR模型的主要区别在于：

无同方差假定：TVP-VAR模型不假定经济变量的同方差结构，更符合实际经济数据的特点

随机波动率捕捉时变影响：模型引入随机波动率，能够更好地捕捉经济变量在不同时期所具有的关系和特征

估计优势：将随机波动率纳入TVP估算中能够显著提高估计性能

在TVP-VAR模型中，所有参数遵循一阶游走过程。随机波动率的概念在TVP-VAR模型中尤为重要，其在经济计量学中发展起来并逐渐应用于宏观经济分析。许多经济数据的产生过程中都伴随着漂移系数和随机波动冲击。为了避免因忽略波动性变化带来的估计偏差，TVP-VAR模型假设了随机波动率。尽管此类模型的似然函数处理复杂，但可以通过贝叶斯推断中的马尔可夫链蒙特卡罗（MCMC）方法进行估计。

二、TVP-VAR建模步骤与注意事项

1. 数据平稳性检验

在进行TVP-VAR建模之前，需确保数据的平稳性。可以采用ADF单位根检验法进行检验。如果数据不平稳，需进行差分处理，直至数据平稳后再进行建模。

2. 变量顺序的影响

变量之间的顺序会影响实证结果的表现。建议将关注的变量放在首位，便于时变关系图的展示和分析。

3. MCMC算法概述

TVP-VAR模型的参数估计通常采用马尔可夫链蒙特卡罗算法（MCMC）。具体步骤包括：

先验密度：π(θ)定义为先验分布

后验分布：π(θ,α,h|y)⁴，涉及时变参数和随机波动率的共轭后验

三、TVP-VAR模型源代码实现

以下为TVP-VAR模型的一个实现示例：

clear; clc;
% 加载数据
load ydata.dat;
load yearlab.dat;
% 模型参数
tau = 40; % 训练样本大小
p = 2; % VAR模型的滞后量
% 生成Z矩阵
ylag = mlag2(Y, p);
ylag = ylag(p + tau + 1:t, :);
K = M + p * (M^2); % 状态向量元素数量
% 初始化模型参数
B0_prmean = B_OLS; B0_prvar = 4 * VB_OLS;
Q_prmean = (0.01)^2 * tau * VB_OLS; Q_prvar = tau;
Sigma_prmean = eye(M); Sigma_prvar = M + 1;
% MCMC参数
nrep = 5000; nburn = 2000; it_print = 100;
%% Gibbs采样
for irep = 1:nrep + nburn
    if mod(irep, it_print) == 0
        disp(irep); toc;
    end
    % 1. 从先验分布采样B_t
    [Btdraw] = carter_kohn_hom(y, Z, Sigmadraw, Qdraw, K, M, t, B0_prmean, B0_prvar);
    % 2. 从逆 Wishart 分布采样Q
    Btemp = Btdraw(:, 2:t)' - Btdraw(:, 1:t-1)';
    sse_2Q = sum(sse_2Q + Btemp(i,:)' * Btemp(i,:), i = 1:t-1);
    Qinv = inv(sse_2Q + Q_prmean);
    Qdraw = wish(Qinv, t + Q_prvar);
    % 3. 从逆 Wishart 分布采样Sigma
    resids = zeros(M, t);
    for i = 1:t
        resids(:, i) = y(:, i) - Z((i-1)*M + 1:i*M, :) * Btdraw(:, i);
    end
    sse_2S = sum(sse_2S + resids(:, i)*resids(:, i)', i = 1:t);
    Sigmainv = inv(sse_2S + Sigma_prmean);
    Sigmadraw = inv(wish(Sigmainv, t + Sigma_prvar));
    % IRF计算
    if irep > nburn
        Bt_postmean = Bt_postmean + Btdraw;
        biga = zeros(M*p, M*p);
        for j = 1:p-1
            biga(j*M + 1:M*(j+1), (j-1)*M + 1:j*M) = eye(M);
        end
        for i = 1:t
            bbtemp = Btdraw(M + 1:K, i);
            for ii = 1:p
                for iii = 1:M
                    biga(iii, (ii-1)*M + 1:ii*M) = bbtemp(splace + 1:splace + M, 1)';
                end
                splace += M;
            end
        end
        shock = eye(M);
        impresp = zeros(M, M*nhor);
        impresp(1:M, 1:M) = shock;
        for j = 1:nhor-1
            impresp(:, j*M + 1:(j+1)*M) = bigj * bigai * bigj' * shock;
            bigai = bigai * biga;
        end
        if yearlab(i, 1) == 1975.00
            impf_m = zeros(M, nhor);
            jj = 0;
            for ij = 1:nhor
                jj += M;
                impf_m(:, ij) = impresp(:, jj);
            end
            imp75(irep - nburn, :, :) = impf_m;
        end
        if yearlab(i, 1) == 1981.50
            impf_m = zeros(M, nhor);
            jj = 0;
            for ij = 1:nhor
                jj += M;
                impf_m(:, ij) = impresp(:, jj);
            end
            imp81(irep - nburn, :, :) = impf_m;
        end
        if yearlab(i, 1) == 1996.00
            impf_m = zeros(M, nhor);
            jj = 0;
            for ij = 1:nhor
                jj += M;
                impf_m(:, ij) = impresp(:, jj);
            end
            imp96(irep - nburn, :, :) = impf_m;
        end
    end
end
% 后续结果分析与可视化
Bt_range = ones(21, 1);
for i = 1:21
    Bt_range(i) = (max(Bt_postmean(i, :)) - min(Bt_postmean(i, :))) / mean(Bt_postmean(i, :));
end
Bt_range = reshape(Bt_range, 3, 7);
rowNames = {'Inflation', 'Unemployment', 'Interest Rate'};
colNames = {'Intercept', 'Inf_1', 'Unemp_1', 'IR_1', 'Inf_2', 'Unemp_2', 'IR_2'};
pc_change_table = array2table(Bt_range, 'RowNames', rowNames, 'VariableNames', colNames);
writetable(pc_change_table, 'pc_change.csv');
% IRF可视化
qus = [0.16, 0.5, 0.84];
imp75XY = squeeze(quantile(imp75, qus));
imp81XY = squeeze(quantile(imp81, qus));
imp96XY = squeeze(quantile(imp96, qus));
figure
set(0, 'DefaultAxesColorOrder', [0 0 0], 'DefaultAxesLineStyleOrder', '--|-|--');
subplot(3, 3, 1)
plot(1:nhor, squeeze(imp75XY(:, 1, :)));
title('Impulse Response of Inflation, 1975:Q1');
ylim([-0.2, 0.1]);
set(gca, 'XTick', 0:3:nhor);
subplot(3, 3, 2)
plot(1:nhor, squeeze(imp75XY(:, 2, :)));
title('Impulse Response of Unemployment, 1975:Q1');
ylim([-0.2, 0.2]);
set(gca, 'XTick', 0:3:nhor);
subplot(3, 3, 3)
plot(1:nhor, squeeze(imp75XY(:, 3, :)));
title('Impulse Response of Interest Rate, 1975:Q1');
ylim([-0.3, 0.1]);
set(gca, 'XTick', 0:3:nhor);
subplot(3, 3, 4)
plot(1:nhor, squeeze(imp81XY(:, 1, :)));
title('Impulse Response of Inflation, 1981:Q3');
ylim([-0.2, 0.1]);
set(gca, 'XTick', 0:3:nhor);
subplot(3, 3, 5)
plot(1:nhor, squeeze(imp81XY(:, 2, :)));
title('Impulse Response of Unemployment, 1981:Q3');
ylim([-0.2, 0.2]);
set(gca, 'XTick', 0:3:nhor);
subplot(3, 3, 6)
plot(1:nhor, squeeze(imp81XY(:, 3, :)));
title('Impulse Response of Interest Rate, 1981:Q3');
ylim([-0.4, 0.1]);
set(gca, 'XTick', 0:3:nhor);
subplot(3, 3, 7)
plot(1:nhor, squeeze(imp96XY(:, 1, :)));
title('Impulse Response of Inflation, 1996:Q1');
ylim([-0.2, 0.1]);
set(gca, 'XTick', 0:3:nhor);
subplot(3, 3, 8)
plot(1:nhor, squeeze(imp96XY(:, 2, :)));
title('Impulse Response of Unemployment, 1996:Q1');
ylim([-0.2, 0.2]);
set(gca, 'XTick', 0:3:nhor);
subplot(3, 3, 9)
plot(1:nhor, squeeze(imp96XY(:, 3, :)));
title('Impulse Response of Interest Rate, 1996:Q1');
ylim([-0.3, 0.1]);
set(gca, 'XTick', 0:3:nhor);
disp('To plot impulse responses, use: plot(1:nhor, squeeze(impXXXY(:, VAR, :))) where VAR=1 for inflation, 2 for unemployment, and 3 for interest rate');

四、备注

完整代码或代写请添加至[QQ2449341593]。

转载地址：http://xlmf.baihongyu.com/

你可能感兴趣的文章