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    endend% 后续结果分析与可视化Bt_range = ones(21, 1);for i = 1:21    Bt_range(i) = (max(Bt_postmean(i, :)) - min(Bt_postmean(i, :))) / mean(Bt_postmean(i, :));endBt_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));figureset(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/

你可能感兴趣的文章