1 时变参数向量自回归（TVP-VAR）介绍
普通向量自回归模型可以研究两个变量之间的相互作用关系，采用脉冲响应可以绘制出某一个变量受到冲击时，另一个变量的响应路径。对一组时间序列数据来说，普通向量自回归只能给出一个脉冲响应函数，但现实中，时间序列数据往往会发生结构性变化，如2008年的金融危机，2020年的新冠疫情，不同时间点上变量之间的脉冲响应路径可能是不一样的。而时变参数向量自回归模型则可以计算出时间序列数据中所有时点上的脉冲响应函数，进而通过对比不同时点上的脉冲响应路径来判断结构变化的程度。
2 问题描述
这里收集了2001年1月到2022年9月的cpi和m2同比增长率数据，采用TVP-VAR模型估计cpi和m2之间的关系，并给出2008年金融危机和2020年新冠疫情后cpi对m2的脉冲响应图。
3 采用R来估计TVP-VAR模型
3.1 安装bvarsv包
bvarsv包中包含估计tvp-var模型的函数。
install.packages("bvarsv")
library(bvarsv)
3.2 导入数据
其中Excel格式数据如下：
安装readxl包，并将excel中的数据导入RStudio
install.packages("readxl")
library(readxl)
mdata=read_excel("D:/M2和CPI.xlsx",sheet="Sheet1")
3.3 设置m2和cpi变量
tsm2=ts(mdata$M2)
tscpi=ts(mdata$CPI)
tms=cbind(tsm2,tscpi)
3.4 估计TVP-VAR模型
bvar.sv.tvp()函数用于估计TVP-VAR模型，其中参数Y是建模的时间序列变量，参数p是滞后期，参数tau为训练集，在脉冲响应分析时需要将其去掉，参数nrep是MCMC重复次数，默认为50000次，参数nburn为燃烧次数，默认为5000次。
本例中设置滞后期为6期，训练集为前24期；且仅为了演示，所以nrep和nburn参数设定的次数比较少。
result=bvar.sv.tvp(Y=tms[,c(1,2)],p=6,tau=24,nrep=5000,nburn=500)
3.5 估计2008年9月cpi对m2的脉冲响应
2008年9月发生了全球金融危机，所以估计该期的脉冲响应。
impulse.responses()函数用于计算脉冲响应，其中参数result为bvar.sv.tvp()函数的返回值，1和2分别为冲击变量和响应变量，t为时期，nhor为脉冲响应的滞后期，scenario为估计脉冲响应的方法，可以取值1，2，3三种方法。
本例中，1表示m2，2表示cpi，2008年9月对应原始序列的69个数（注意原始序列中2008年9月是第93个数，需要把前24个数丢弃，所以变为69）；脉冲响应滞后期设为30。
b1=impulse.responses(result,1,2,t=69,nhor=30,scenario=1)
从估计的脉冲响应图可以看出，m2发生冲击后，cpi前5期处于波动状态，第6期至8期开始上升，之后逐渐恢复。
3.6 估计2020年1月cpi对m2的脉冲响应
2020年1月发生了新冠疫情，所以估计该期的脉冲响应。
本例中，2020年1月对应原始序列的205个数（注意原始序列中2020年1月是第229个数，需要把前24个数丢弃，所以变为205）；其余参数设置与上面相同。
b2=impulse.responses(result,1,2,t=205,nhor=30,scenario=1)
从估计的脉冲响应图可以看出，m2发生冲击后，前16期cpi在缓慢上升，之后逐渐恢复。与2008年8月的脉冲响应相比，2020年1月份cpi对m2的脉冲响应持续期更久。
3.7 获取脉冲响应的具体数值
impulse.responses()函数的返回值中包含了500组脉冲响应的序列，上面两幅图中的脉冲响应路径是这500组脉冲响应的中位数组成的序列。有些时候我们需要获取该脉冲响应序列。
b1_irf和b2_irf即为获取的脉冲响应值。
install.packages("robustbase")
library("robustbase")
b1_irf=colMedians(b1$irf)
b2_irf=colMedians(b2$irf)
将b1_irf和b2_irf脉冲响应绘制在一副图上。
plot(b2_irf,type="o",main = "cpi对m2的脉冲响应")
lines(b1_irf,type="o",col="blue")