README.md
In juhokalle/rmfd4dfm: Estimation of DFMs using RMFD-E parametrization

Estimation of Impulse-Response Functions with Dynamic Factor Models: A New Parametrization ================

This repository contains the associated R-package and data to the paper Estimation of Impulse-Response Functions with Dynamic Factor Models: A New Parametrization authored by Juho Koistinen and Bernd Funovits (2022). Specifically, the functions within the package allow the user to estimate a structural dynamic factor model (DFM), with the common component parametrized as a right matrix fraction description in echelon form (RMFD-E), which is introduced in Section 2.2 of the associated paper. This document also shows how to replicate the empirical exercise of Section 4.

The package builds on the R-packages rationalmatrices and RLDM, authored by Bernd Funovits and Wolfgang Scherrer. Since these packages might change, the parts which are necessary for the analysis in the associated article are extracted to R files \~/rmfd4dfm/R/zz_ratmat and \~/rmfd4dfm/R/zz_rldm. Importantly, please reference the rationalmatrices and the RLDM packages should you use their functionalities and not this package.

The package can be installed using the following command as is the standard with packages hosted by GitHub.

# install.packages("devtools")
devtools::install_github("juhokalle/rmfd4dfm")

Here I will provide a brief and hopefully accessible introduction to the modelling approach presented in the paper, while the details are available from the paper. The aim is to identify and estimate impulse-response functions (IRFs) using dynamic factor models (DFMs), with the emphasis placed on macroeconomic applications and analysis. The model of interest is given as

$(1): \quad x_{t} =d_{0}z_{t}^{*}+d_{1}z_{t-1}^{*}+\cdots+d_{s}z_{t-s}^{*}+\xi_{t}$

$(2): \quad z_{t}^{*} =c_{1}z_{t-1}^{*}+\cdots+c_{p}z_{t-p}^{*}+\varepsilon_{t}$

$(3): \quad \varepsilon_{t}=Hu_{t}$

where

$x_t$ is an $n$ -dimensional observed time series, with usually $n>100$ .
$z_t^*$ is a $q$ -dimensional dynamic factor process, with usually $q\<10$ .
$x_t$ is an $n$ -dimensional idiosyncratic term, we assume $\mathbb E(\xi_t \xi_t')=\sigma_\xi^2I_n$ .
$\varepsilon_t$ is the innovation to $z_t^*$ , with $\mathbb E(\varepsilon_t \varepsilon_t')=\Sigma_\varepsilon$ .
$\{d_i\}_{i=0}^s$ and $\{c_j\}_{i=1}^p$ are $n \times q$ and $q \times q$ parameter matrices, values of which are at the center of interest as the IRFs are constructed from these matrices.
$u_t$ is the $q$ -dimensional structural shock process, with $\mathbb E (u_t u_t') = I_q$ .
$H$ is $(q\times q)$ -dimensional structural impact multiplier matrix, which we identify as $H=chol(\Sigma_\varepsilon)$ , where $chol(\Sigma_\varepsilon)$ is the lower triangular Cholesky factor of $\Sigma_\varepsilon$ . Note that any other identification method for uncovering $H$ is valid as in the structural vector autoregression (SVAR) analysis.

We can write Eq. (2) compactly as $z_t^*=c(L)^{-1}\varepsilon_t$ and substitute this and Eq. (3) into Eq. (1) to get

$x_t=d(L)c(L)^{-1}H u_t + \xi_t,$

where

$c(L)=I_q-c_1 L - \cdots - c_p L^p$

and

$d(L)=d_0 + d_1 L + \cdots + d_sL^s.$

If we assume that the idiosyncratic component accounts for measurements errors or sectoral dynamics that pertain to a small number of variables, we can use the structural IRF $k(L)H$ to study the shock propagation from $q$ structural shocks to $n$ observed macroeconomic variables, with $k(L)=d(L)c(L)^{-1}$ . The IRF $k(L)$ is not identified without further restrictions since one can always post-multiply $d(L)$ and $c(L)$ by some $q\times q$ polynomial matrix $m(L)$ such that it “cancels out”:

$k(L)=d(L)c(L)^{-1}=[d(L)m(L)][c(L)m(L)]^{-1}=\bar d(L) \bar c(L)^{-1}.$

The problem is that the researcher cannot distinguish between $d(L)c(L)^{-1}$ and $\bar d(L) \bar c(L)^{-1}$ from the first and second moments of the data, and therefore attempts at drawing conclusions from the structural IRF $k(z)H$ are futile without further assumptions on $c(L)$ and $d(L)$ .

Our insight is that $k(L)$ can be identified, that is, the set of matrices $m(L)$ can be narrowed down to $I_q$ , using the identification restrictions that are standard in the literature dealing with the identification of the vector autoregressive moving average (VARMA) models (for a recent summary, see Deistler and Scherrer, 2019). In this model class, the IRF is given as $k(L)=a(L)^{-1}b(L)$ for some AR and MA lag polynomials $a(L)$ and $b(L)$ of suitable dimensions, and here $a(L)$ and $b(L)$ can be pre-multiplied by some lag polynomial $m(L)$ to obtain an observationally equivalent IRF:

$k(L)=a(L)^{-1}b(L)=[m(L)a(L)]^{-1}[m(L)b(L)]=\bar a(L)^{-1}\bar b(L).$

One popular identification approach for the VARMA model is to use the canonical echelon form parametrization, which involves zero and unity restrictions in $a(L)$ and $b(L)$ . By noting that the identification restrictions for the IRF in right matrix fraction description (RMFD, for short) model, i.e. $k(L)=d(L)c(L)^{-1}$ , are equivalent to those for the VARMA model after transposing, the derivation is quite straightforward using the existing results for the VARMA model (for details, see Appendix A).

The structure of the RMFD model in echelon form is quite complicated, and while the functions of this package rmfd4dfm code these restrictions without the need for the user to trouble herself with the implementation, one might still wonder what is the point of going through the painstaking process of deriving and coding the restrictions in the first place, as there are alternatives in the literature? To make the point, let us introduce a popular alternative used in the literature and compare it with our approach. One way to deal with the identification and estimation of the DFM presented above is to stack the dynamic factors into a $r=q(s+1)$ -dimensional vector

$z_t = \left(z_t^{*'}, z_{t-1}^{*'}, \cdots, z_{t-s}^{*'} \right)',$

and write Eq. (1) as

$x_t = \left(d_0,\cdots,d_s\right)z_t + \xi_t = Dz_t + \xi_t,$

which makes the model amenable to the principal components (PC) methods. In the second step, the static factor process is modelled as a VAR( $k$ ) process

$z_t - A_1 z_{t-1} - \cdots - A_k z_{t-k} = A(L)z_t = B\varepsilon_t,$

where $A(L)$ is an $r \times r$ VAR lag polynomial of order $k$ , and $B$ is an $r\times q$ constant matrix. After the estimation of $D$ and $A(L)$ , the IRF $k(L)=DA(L)^{-1}B$ can be constructed straightforwardly.

This alternative using PCs is straightforward from an estimation point of view, but let us point two associated restrictive features. First, note that the minimal number of identifying restrictions needed is $r^2=q^2(s+1)^2$ , which is higher than that of in our parametrization, i.e. $q^2$ , whenever $s>0$ (Bai and Wang, 2012). Second, we note that the reliable estimation of VAR on $z_t$ can be difficult if the lag order is misspecified, which is due to the singularity of the covariance matrix of the innovation to $z_t$ , i.e. $\mathbb E (B\varepsilon_t \varepsilon_t'B')=B\Sigma_\varepsilon B'$ has rank $q$ , which is smaller than $r=(s+1)q$ if $s>0$ (Hörmann and Nisol, 2021). On the other hand, the law of motion for $z_t^*$ is standard non-singular VAR in our setup. These two examples make the case that the “static” method is suboptimal whenever $s>0$ and warrant considering other alternatives, development of which is the main goal of our paper.

The package makes available three different macroeconomic panels of data measured at a monthly frequency. First is the data used by Forni and Gambetti (2010), and the documentation can be accessed using command ?FG_data in the R console. Second and third data sets, FRED_light and FRED_heavy originate from the FRED-MD data set, which is documented carefully in McCracken and Ng (2015). The difference between the data sets concerns the transformations to obtain stationarity, while the time span is same in both data sets to facilitate comparison to Forni and Gambetti (2010). The details can be accessed by ?FRED_light or ?FRED_heavy.

The empirical example compares the DFMs estimated by Forni and Gambetti (2010) to that developed in our paper. The idea is to assess the effects of monetary policy to key macroeconomic variables via structural DFMs, which guards against the omitted variable bias, to which the SVAR methods can be more susceptible.

The following code snippets show how to replicate the results given in the empirical section of our paper. First, we order the variables of interest in the structural analysis variables first in the data matrix, i.e. industrial production (INDPRO, the mnemonic in McCracken and Ng, 2015), consumer price index (CPIAUCSL), federal funds rate (FEDFUNDS), and real Swiss/US exchange rate (EXSZUSx), in this order. In the structural analysis with recursive identification, the ordering is important, as it determines the short-run restrictions imposed on the structural impact multiplier matrix. (For a brief and to the point summary of the structural identification using DFMs, see Section 9 of the lecture notes written by Matteo Barigozzi.) Additionally, one must include a vector int_ix in the data object FRED_heavy coding the positions of the variables of interest, which in this case is simply a sequence 1:4.

library("rmfd4dfm")
# code the positions of the variables of interest
int_vars_fred <- c("INDPRO", "CPIAUCSL", "FEDFUNDS", "EXSZUSx")
int_ix <- sapply(int_vars_fred, function(x) which(names(FRED_heavy$df)==x))
# re-organize data matrix s.t. the variables of interest are ordered first
perm_ix <- c(int_ix, (1:ncol(FRED_heavy$df))[-int_ix])
FRED_heavy$df <- FRED_heavy$df[,perm_ix]
FRED_heavy$trans_ix <- FRED_heavy$trans_ix[perm_ix]
FRED_heavy$int_ix <- 1:4

Second, we need to determine the static factor dimension. This is needed for an estimate of the lag structure in eqs. (1)–(2), while in the “static” method it defines the dimension of $z_t$ . To this end, we use functions baingcriterion and abc_crit, which implement the tests developed in Bai and Ng (2002) and Alessi, Barigozzi and Capasso (2009), respectively.

# Estimate the number of static factors ####
baing_cr <- baingcriterion(FRED_heavy$df, rmax = 25)$IC[1:2]
abc_cr <- replicate(50, abc_crit(FRED_heavy$df, kmax = 25))
get_mode <- function(x) unique(x)[which.max(colSums(sapply(unique(x), function(z) x %in% z)))]
r_hats <- c("ICp1" = baing_cr[1],
            "ICp2" = baing_cr[2],
            "ABCp1" = get_mode(unlist(abc_cr[1,])),
            "ABCp2" = get_mode(unlist(abc_cr[2,]))
            )

Upon fixing the state dimension of the state space representation of eqs. (1)–(2), estimation of the DFM can be carried out. Here we fix the static factor dimension to 8. In the RMFD-E parametrization of the DFM and for a given state dimension, it is likely that one has many model specifications available from which to choose “the best model”. Therefore, the function do_everything_rmfd estimates all the model alternatives that are consistent with the model selection criteria given in Section 3.3. and chooses the one that minimizes Bayesian information criterion. For this model, the function then estimates confidence intervals for the IRF using block bootstrap. Finally, the function returns the impulse responses of the variables to the third shock, which is standardized to have an immediate impact of 0.5. For responses to other shocks and/or size, one can include an index vector of length two, which specifies the shock of interest and normalization constant, in the data object specified as the first argument of do_everything_rmfd by df$shock_ix. For a more flexible estimation setup, the user can consult function estim_wrap, which produces an estimate of the non-structural IRF $k(L)=d(L)c(L)^{-1}$ for a given model structure. The functions do_everything_fglr and do_everything_svar perform the same task for the static DFM and structural VAR, respectively.

# D-DFM
est_obj <- do_everything_rmfd(df = FRED_heavy,
                              r = 8,
                              h = 50,
                              nrep = 500,
                              conv_crit = 1e-5,
                              ci = 0.68,
                              init_rep = 1,
                              verbose = TRUE)
# S-DFM
est_fglr <- do_everything_fglr(df = FRED_heavy,
                               r = 8,
                               k = 2,
                               h = 50,
                               nrep = 500,
                               ci = 0.68)

# SVAR
svar_irf <- do_everything_svar(df = FRED_heavy,
                               p = 9,
                               nrep = 500,
                               h = 50,
                               ci = 0.68)

Finally, the plot_irfs function returns the IRFs along with the possible error bands in a nice plot. To gather the IRF plots into one figure, which can then be exported to pdf, for example, we use gridExtra::marrangeGrob.

p1 <- plot_irfs(est_obj$irf, int_vars_fred, "D-DFM")
p2 <- plot_irfs(est_fglr, int_vars_fred, "S-DFM", label_y = FALSE)
p3 <- plot_irfs(svar_irf, int_vars_fred, "SVAR", label_y = FALSE)
plot1 <- gridExtra::marrangeGrob(c(p1,p2,p3), nrow = 4, ncol = 3, top = NULL)

juhokalle/rmfd4dfm documentation built on July 18, 2024, 10:19 p.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

Tweet to @rdrrHQ

GitHub issue tracker

ian@mutexlabs.com