In this example we will model the long run behaviour of nonstationary time series with a Bayesian vector autoregressive model that makes use of a Minnesota prior. The model for this example is contained in the file T9_minn.R
. To start off we can clear all the variables from the current environment and close all the plots.
rm(list = ls())
graphics.off()
Thereafter, we will need to make use of the BMR
package, which has a number of dependencies:
install.packages("Rcpp")
install.packages("ggplot2")
devtools::install_github("kthohr/BMR")
To make use of the routines in the BMR
package we need to make use of the library
command.
library(BMR)
As this data is contained in a .csv
file we need to set the directory to tell R where to find the datafile. You will need to change the following extension to ensure that the correct path is specified.
setwd("C:\\Users\\image")
This allows us to load the data and transform the price index to yearonyear inflation, we proceed as follows.
dat < read.csv(file = "USdata.csv")
To ensure that this has all been completed correctly, we can then plot the data, using the gtsplot
command from within BMR
.
gtsplot(dat[, 2:4], dates = dat[, 1])
This graph would suggest that all the variables would potentially have a stochastic trend.
We firstly need to create a model object that will store all the results from the estimation procedure. We can call this object bvar_obj
.
bvar_obj < new(bvarm)
The basic structure of the model and the data need to be included in the build
part of this object. In this case we are going to make use of a VAR(4) with a constant and the data that we will be using is stored in columns two through four.
bvar_obj$build(data.matrix(dat[,2:4]),
TRUE, # constant
4) # lags
To use randomwalk priors for the nonstationary variables, we set the prior mean values to unity before estimating the model parameters.
prior < c(1, 1, 1)
To complete the construction of the prior we need to select values for the hyperparameters. In this case we are going to follow @Canova:2007 and set the first hyperparameter, \(\lambda\), to 0.2 and the second hyperparameter, \(\theta\), to 0.5. He also recommends the values that were used for the second and third hyperparameters, HP2
and HP3
. The other elements take on the default vaules that are taken from @Koop:2010.
bvar_obj$prior(prior, # prior mean value
1, # var_type
1, # decay_type
0.2, # HP1
0.5, # HP2
10^5, # HP3
1.0) # HP4
The final requirement is to specify the number of draws that are required for the Gibbs sampler, which we are setting to 10,000
in this case.
bvar_obj$gibbs(10000)
To plot the density functions for the parameter estimates we would then make use of the following commands.
plot(bvar_obj, var_names = colnames(dat)[1], save = FALSE)
where the mean and variances for this distributions are stored in the objects, bvar_obj$alpha_pt_mean
and bvar_obj$alpha_pt_var
. And then the impulse response functions could be plotted as follows:
IRF(bvar_obj, 20, var_names = colnames(dat)[1], save = FALSE)
To generate forecasts for each of the variables then we proceed as follows, where shocks
indicates that we would like to incorporate the effect of the shocks in estimation. The amount of back_data
refers to the amount of observations to include in the graph prior to the initial forecasting point.
forecast(bvar_obj, shocks = TRUE, var_names = colnames(dat)[1],
back_data = 12, save = FALSE)
Note that the forecasts are future estimates of the trend in the process, and we are no longer forecasting time taken to return to trend. The specific values for the forecasts of the individual variables can be extracted from the object by amending the object to:
predict < forecast(bvar_obj, shocks = TRUE, var_names = colnames(dat)[1],
save = TRUE)
The model for this example is contained in the file T9_sign.R
. To start off we can clear all the variables from the current environment and close all the plots.
rm(list = ls())
graphics.off()
Thereafter, we will need to make use of the VARsignR
package that has been installed. If you need install this package, whcih is available on cran then you would want to run the following routine: install.packages("VARsignR")
.
library(VARsignR)
As this data is contained in a .csv
file we need to set the directory to tell R where to find the datafile. You will need to change the following extension to ensure that the correct path is specified.
setwd("C:\\Users\\image")
This allows us to load the data and transform the price index to yearonyear inflation, we proceed as follows.
dat_raw < read.csv(file = "uhlig_data.csv")
dat < as.ts(dat_raw[, 2:7], start = c(1959, 1), frequency = 12)
To ensure that this has all been completed correctly, we can then plot the data, using the gtsplot
command from within BMR
.
plot(dat)
In terms of the restrictions we want to impose the following:
Hence the constraint would take the form:
constr < c(+4, 3, 2, 5)
To replicate the results of @Uhlig:2005 we would then make use of a VAR(12) and no constant in the model. To generate the impulse response functions we make use of 60 steps. We use 200 draws from the posterior and 200 subdraws for each posterior draw to generate the impulse vectors and the candidate impulse responses to which the rejection algorithm will be applied.
model1 < uhlig.reject(Y = dat, nlags = 12, draws = 200,
subdraws = 200, nkeep = 1000, KMIN = 1, KMAX = 6, constrained = constr,
constant = FALSE, steps = 60)
##

  0%

  1%

=  2%

=  3%

==  4%

==  5%

==  6%

===  6%

===  7%

===  8%

====  8%

====  9%

====  10%

=====  10%

=====  11%

=====  12%

======  12%

======  13%

======  14%

=======  14%

=======  15%

=======  16%

========  17%

========  18%

=========  19%

=========  20%

=========  21%

==========  22%

==========  23%

===========  24%

===========  25%

===========  26%

============  26%

============  27%

============  28%

=============  28%

=============  29%

=============  30%

==============  30%

==============  31%

==============  32%

===============  32%

===============  33%

===============  34%

================  34%

================  35%

================  36%

=================  37%

=================  38%

==================  39%

==================  40%

==================  41%

===================  42%

===================  43%

====================  44%

====================  45%

====================  46%
summary(model1)
## Length Class Mode
## IRFS 360000 none numeric
## FEVDS 360000 none numeric
## SHOCKS 456000 none numeric
## BDraws 38880 none numeric
## SDraws 3240 none numeric
The results of the model are then provided as a result of the summary command, which gives an indication of what has been stored in the object model1
. To plot the impulse response functions we would first give the variables names before creating an object irfs1
. Thereafter, we would make use of the irfplot
command to display the result:
nam < c("GDP", "GDP Deflator", "Comm Price Index", "Fed Funds Rate",
"NB Reserves", "Total Reserves")
irfs1 < model1$IRFS
irfplot(irfdraws = irfs1, type = "median", labels = nam,
save = FALSE, bands = c(0.16, 0.84), grid = TRUE, bw = FALSE)
Note that in this case we have used confidence intervals of 16%. We could also plot the forecasterror variance decomposition using a similar procedure, where the object fevd1
contains the results and fevdplot
displays these results in a graph.
fevd1 < model1$FEVDS
fevdplot(fevd1, label = nam, save = FALSE, bands = c(0.16,
0.84), grid = TRUE, bw = FALSE, table = FALSE, periods = NULL)
As an alternative we could display the results in a table, by making use of the following commands:
fevd.table < fevdplot(fevd1, table = TRUE, label = nam,
periods = c(1, 10, 20, 30, 40, 50, 60))
print(fevd.table)
## GDP GDP Deflator Comm Price Index Fed Funds Rate
## 1 8.86 6.91 7.18 10.17
## 10 9.52 8.16 8.00 12.53
## 20 12.06 9.22 8.33 12.96
## 30 13.11 9.46 8.72 12.93
## 40 12.67 9.62 9.22 12.84
## 50 12.40 9.79 9.56 12.72
## 60 12.19 9.94 9.93 12.66
## NB Reserves Total Reserves
## 1 5.72 7.73
## 10 8.82 8.57
## 20 10.23 9.57
## 30 10.26 9.75
## 40 10.34 9.95
## 50 10.49 10.02
## 60 10.57 10.17
To recover the shocks we would need to extract the relevant information that is stored in the model object. Note that we have kept 1000 simulations, so we use the quantile
function to generate confidence intervals for these shocks. These may be stored as a time series object, where we lost the first 12 observations as a part of the estimation procedure.
shocks < model1$SHOCKS
ss < ts(t(apply(shocks, 2, quantile, probs = c(0.5, 0.16,
0.84))), frequency = 12, start = c(1960, 1))
To plot the result, we could make use of the following commands:
plot(ss[, 1], type = "l", col = "blue", ylab = "Interest rate shock",
ylim = c(min(ss), max(ss)))
abline(h = 0, col = "black")
lines(ss[, 2], col = "red")
lines(ss[, 3], col = "red")
model2 < rwz.reject(Y = dat, nlags = 12, draws = 200, subdraws = 200,
nkeep = 1000, KMIN = 1, KMAX = 6, constrained = constr,
constant = FALSE, steps = 60)
##

  0%

  1%

=  2%

=  3%

==  4%

==  5%

==  6%

===  6%

===  7%

===  8%

====  8%

====  9%

====  10%

=====  10%

=====  11%

=====  12%

======  12%

======  13%

======  14%

=======  14%

=======  15%

=======  16%

========  17%

========  18%

=========  19%

=========  20%

=========  21%

==========  22%

==========  23%

===========  24%

===========  25%

===========  26%

============  26%

============  27%

============  28%

=============  28%

=============  29%

=============  30%

==============  30%

==============  31%

==============  32%

===============  32%

===============  33%

===============  34%

================  34%

================  35%

================  36%

=================  37%

=================  38%

==================  39%

==================  40%

==================  41%

===================  42%

===================  43%

====================  44%

====================  45%

====================  46%

=====================  46%

=====================  47%

=====================  48%

======================  48%

======================  49%

======================  50%

=======================  50%

=======================  51%

=======================  52%

========================  52%

========================  53%

========================  54%

=========================  54%

=========================  55%

=========================  56%

==========================  57%

==========================  58%

===========================  59%

===========================  60%

===========================  61%

============================  62%

============================  63%

=============================  64%

=============================  65%

=============================  66%

==============================  66%

==============================  67%

==============================  68%

===============================  68%

===============================  69%

===============================  70%

================================  70%

================================  71%

================================  72%

=================================  72%

=================================  73%

=================================  74%

==================================  74%

==================================  75%

==================================  76%

===================================  77%

===================================  78%

====================================  79%

====================================  80%

====================================  81%

=====================================  82%

=====================================  83%

======================================  84%

======================================  85%

======================================  86%

=======================================  86%

=======================================  87%
irfs2 < model2$IRFS
irfplot(irfdraws = irfs2, type = "median", labels = nam,
save = FALSE, bands = c(0.16, 0.84), grid = TRUE, bw = FALSE)
model3 < uhlig.penalty(Y = dat, nlags = 12, draws = 2000,
subdraws = 1000, nkeep = 1000, KMIN = 1, KMAX = 6, constrained = constr,
constant = FALSE, steps = 60, penalty = 100, crit = 0.001)
##

  0%

  1%

=  1%

=  2%

=  3%

==  3%

==  4%

==  5%

==  6%

===  6%

===  7%

===  8%

====  8%

====  9%

====  10%

=====  10%

=====  11%

=====  12%

======  12%

======  13%

======  14%

=======  14%

=======  15%

=======  16%

=======  17%

========  17%

========  18%

========  19%

=========  19%

=========  20%

=========  21%

==========  21%

==========  22%

==========  23%

===========  23%

===========  24%

===========  25%

===========  26%

============  26%

============  27%

============  28%

=============  28%

=============  29%

=============  30%

==============  30%

==============  31%

==============  32%

===============  32%

===============  33%

===============  34%

================  34%

================  35%

================  36%

================  37%

=================  37%

=================  38%

=================  39%

==================  39%

==================  40%

==================  41%

===================  41%

===================  42%

===================  43%

====================  43%

====================  44%

====================  45%

====================  46%

=====================  46%

=====================  47%

=====================  48%

======================  48%

======================  49%

======================  50%
irfs3 < model3$IRFS
irfplot(irfdraws = irfs3, type = "median", labels = nam,
save = FALSE, bands = c(0.16, 0.84), grid = TRUE, bw = FALSE)
constr5 < c(+4, 3, 2, 5, 6)
model5c < uhlig.reject(Y = dat, nlags = 12, draws = 200,
subdraws = 200, nkeep = 1000, KMIN = 1, KMAX = 6, constrained = constr5,
constant = FALSE, steps = 60)
##

  0%

  1%

=  2%

=  3%

==  4%

==  5%

==  6%

===  6%

===  7%

===  8%

====  8%

====  9%

====  10%

=====  10%

=====  11%

=====  12%

======  12%

======  13%

======  14%

=======  14%

=======  15%

=======  16%

========  17%

========  18%

=========  19%

=========  20%

=========  21%

==========  22%

==========  23%

===========  24%

===========  25%

===========  26%

============  26%

============  27%

============  28%

=============  28%

=============  29%

=============  30%

==============  30%

==============  31%

==============  32%

===============  32%

===============  33%

===============  34%

================  34%

================  35%

================  36%

=================  37%

=================  38%

==================  39%

==================  40%

==================  41%

===================  42%

===================  43%

====================  44%

====================  45%

====================  46%

=====================  46%

=====================  47%

=====================  48%

======================  48%

======================  49%

======================  50%

=======================  50%

=======================  51%

=======================  52%

========================  52%

========================  53%

========================  54%

=========================  54%

=========================  55%

=========================  56%

==========================  57%

==========================  58%

===========================  59%

===========================  60%

===========================  61%

============================  62%

============================  63%

=============================  64%
irf5c < model5c$IRFS
irfplot(irf5c, labels = nam)
fp.target(Y = dat, irfdraws = irfs1, nlags = 12, constant = F,
labels = nam, target = TRUE, type = "median", bands = c(0.16,
0.84), save = FALSE, grid = TRUE, bw = FALSE, legend = TRUE,
maxit = 1000)