MultiATSM: An R Package for Arbitrage-Free Macrofinance Multicountry Affine Term Structure Models

The MultiATSM package provides estimation tools and a wide range of outputs for eight macrofinance affine term structure model (ATSM) classes, supporting practitioners, academics, and policymakers. All models extend the single-country framework of Joslin et al. (2014) to multicountry settings, with additional adaptations from Jotikasthira et al. (2015) and Candelon and Moura (2024). These model extensions incorporate, respectively, the presence of a dominant (global) economy and adopt a global vector autoregressive (GVAR) setup to capture the joint dynamics of risk factors. The package generates diverse outputs for each ATSM, including graphical representations of model fit, risk premia, impulse response functions, and forecast error variance decompositions. It also implements bootstrap methods for confidence intervals and produces bond yield forecasts.

1 Introduction

The term structure of interest rates (or yield curve) describes the relationship between bond yields and investment maturities. As Piazzesi (2010) emphasizes, understanding its dynamics is essential for several reasons. First, long-term yields incorporate market expectations of future short-term rates, making the yield curve a handy forecasting tool for macroeconomic aggregates like output and inflation. As such, this supports optimal consumption-saving decisions and capital allocation by economic agents. Second, it plays a key role in the transmission of monetary policy, linking short-term policy rates to long-term borrowing costs. Third, it guides fiscal authorities in shaping debt maturities to balance refinancing risk and interest rate exposure. Fourth, it is essential for pricing and hedging interest rate derivatives, which rely on accurate yield curve modelling.

Affine Term Structure Models (ATSMs) are the workhorse in yield curve modelling. Based on the assumption of no arbitrage, ATSMs offer a flexible framework to assess how investors price risks and generate predictions for the price of any bond (see Piazzesi (2010); Gürkaynak and Wright (2012) for comprehensive reviews). Early ATSMs gained popularity for their ability to capture nearly all term structure fluctuations, appealing to both academics and practitioners (Vasicek 1977; Duffie and Kan 1996; Dai and Singleton 2002). While these models produce accurate statistical descriptions of the yield curve, they are silent on the deeper economic determinants that policymakers require for causal inference.

In response to this limitation, a large body of research has emerged to explore the interplay between the term structure and macroeconomic developments (seminal contributions include Ang and Piazzesi (2003) and Rudebusch and Wu (2008)). A prominent contribution in this area is the unspanned economic risk framework developed by Joslin et al. (2014) (henceforth JPS, 2014). In essence, this model assumes an arbitrage-free bond market and considers a linear state space representation to describe the dynamics of the yield curve. Compared to earlier macrofinance ATSMs, JPS (2014) offers a tractable estimation approach that integrates traditional yield curve factors (spanned factors) with macroeconomic variables (unspanned factors). As a result, the model delivers a strong cross-sectional fit while explicitly linking bond yield responses to the state of the economy.

The work of JPS (2014) lays the foundational framework for the modelling tools included in the MultiATSM package (Moura 2025). In addition to the original single-country setup proposed by JPS (2014), the package incorporates multicountry extensions developed by Jotikasthira et al. (2015) (henceforth JLL, 2015) and Candelon and Moura (2024) (henceforth CM, 2024). Altogether, the package offers functions to build eight types of ATSMs, covering the original versions and several variants of these three frameworks.

Beyond complete routines for model estimation, MultiATSM produces a wide range of analytical outputs. In particular, it generates graphical representations such as model-implied bond yields, bond risk premia, and both orthogonalized and generalized versions of: (i) impulse response functions, and (ii) forecast error variance decompositions for yields and risk factors. Confidence intervals for the two latter outputs can be computed using three bootstrap methods: residual-based, block, or wild bootstrap. Moreover, the package supports out-of-sample forecasting of bond yields across the maturity spectrum. This paper provides detailed guidance on how to use the MultiATSM package effectively.

There are a few notable packages for term structure modelling in the R programming environment. YieldCurve (Guirreri 2015) and fBonds (Setz 2017) provide a collection of functions to build term structures based on the frameworks of Nelson and Siegel (1987) and Svensson (1994). These yield curve methods have gained popularity for their parsimonious parameterization and good empirical fit. However, these models do not rule out arbitrage opportunities, a limitation addressed by ATSMs. Moreover, the focus of YieldCurve and fBonds is restricted to parameter estimation and yield curve fitting, without offering additional model outputs such as those provided by MultiATSM.

Several other R packages support time series modeling (Hyndman and Killick 2025), particularly within state space and vector autoregressive (VAR) frameworks. State space packages are relatively few and tend to focus on either estimation, statespacer (Beijers 2023), or simulation, simStateSpace(Pesigan 2025). VAR-based tools are more numerous. For instance, vars (Pfaff and Stigler 2024) and MTS (Tsay et al. 2022) provide extensive functionality for estimation, diagnostics, and forecasting, while svars (Lange et al. 2023) adds structural identification methods. High-dimensional VARs are handled by packages like bigtime (Wilms et al. 2023) and BigVAR (Nicholson et al. 2025), and cross-country spillovers are modeled by Spillover (Urbina 2024) and BGVAR (Boeck et al. 2024).

Although these tools share some features with MultiATSM, they are tailored to standard state space or VAR analysis. In contrast, MultiATSM embeds VAR dynamics within a state space representation that is explicitly grounded in arbitrage-free asset pricing theory. As such, MultiATSM fills a specific gap in the R ecosystem by combining the structure of ATSMs with the flexibility of modern time series tools.

The remainder of the paper is organized as follows. Section 2 outlines the theoretical foundations of the ATSMs implemented in the MultiATSM package, and Section 3 details each model’s features. The subsequent sections focus on the practical implementation of ATSMs. Section 4 presents the dataset included in the package. Section 5 explains the user inputs required for model estimation. Section 6 explains the estimation procedure, and Section 7 shows how to estimate ATSMs from scratch using MultiATSM. Replications of published academic studies are provided in the Appendix.

2 ATSMs with unspanned economic risks: theoretical background

In this section, I outline several arbitrage-free ATSMs with unspanned macroeconomic risks available in the MultiATSM package. A key appealing feature of these setups is their ability to disentangle the yield curve into a cross-sectional component, governed by the risk-neutral (\(\mathbb{Q}\)) dynamics, and a time-series component, driven by the physical (\(\mathbb{P}\)) dynamics. In light of this characteristic of the models, I present the single and the multicountry \(\mathbb{Q}\)-dynamics model dimensions in Section 2.1. Next, I expose the specific features of the risk factor dynamics under the \(\mathbb{P}\)-measure of the various restricted and unrestricted VARs settings in Section 2.2. Section 2.3 describes the model estimation procedures.

Model cross-sectional dimension (Q-dynamics)

Single-country specifications (individual Q-dynamics model classes)

The model cross-sectional structure is based on two central equations. The first one assumes that the country \(i\) short-term interest rate at time \(t\), \(r_{i,t}\), is an affine function of \(N\) unobserved (latent) country-specific factors, \(\boldsymbol{X_{i,t}}\): \[\begin{equation} \underset{(1 \times 1)}{\vphantom{\Big|} r_{i,t}} = \underset{(1 \times 1)}{ \vphantom{\Big|} \delta_{i,0}} + \underset{(1 \times N)}{% \vphantom{\Big|} \boldsymbol{\delta}_{i,1}^{\top}} \underset{(N \times 1)}{% \vphantom{\Big|} \boldsymbol{X}_{i,t}}\text{,} \tag{1} \end{equation}\] where \(\delta_{i,0}\) and \(\boldsymbol{\delta_{i,1}}\) are time-invariant parameters.

The second equation assumes that the unobserved factor dynamics for each country \(i\) follow a maximally flexible, first-order, \(N-\)dimensional multivariate Gaussian (\(\mathcal{N}\)) VAR model under the \(\mathbb{Q}\)-measure:

\[\begin{align} & \underset{(N \times 1)}{\boldsymbol{\vphantom{\Big|} X_{i,t}}} = \underset{(N \times 1)}{\boldsymbol{\vphantom{\Big|} \mu^{Q}_{i,X}}} + \underset{(N \times N)}{\vphantom{\Big|} \Phi^{Q}_{i,X}} \underset{(N \times 1)}{\boldsymbol{\vphantom{\Big|} X_{i,t-1}}} + \underset{(N \times N)}{\vphantom{\Big|} \Gamma_{i,X}} \underset{(N \times 1)}{\boldsymbol{\vphantom{\Big|} \varepsilon_{i,t}^{Q}}}\text{,} & \boldsymbol{\varepsilon_{i,t}^{Q}}\sim {\mathcal{N}_N}(\boldsymbol{0}_N,\mathrm{I}_N)\text{,} \tag{2} \end{align}\] where \(\boldsymbol{\mu^{Q}_{i,X}}\) contains intercepts, \(\Phi^{Q}_{i,X}\), the feedback matrix, and \(\Gamma_{i,X}\) a lower triangular matrix.

Based on Equations (1) and (2), Dai and Singleton (2000) show that the country-specific zero-coupon bond yield with maturity of \(n\) periods, \(y_{i,t}^{(n)}\), is affine in \(\boldsymbol{X_{i,t}}\): \[\begin{equation} \underset{(1 \times 1)}{\vphantom{\Big|} y_{i,t}^{(n)}} = \underset{(1 \times 1)}{\vphantom{\Big|} a_{i,n}(\Theta_{n})} + \underset{(1 \times N)}{\vphantom{\Big|} \boldsymbol{b_{i,n}(\Theta_{n})}^\top} \underset{(N \times 1)}{\vphantom{\Big|} \boldsymbol{X_{i,t}}}\text{,} \tag{3} \end{equation}\] where \(a_{i,n}(\Theta _{n})\) and \(\boldsymbol{b_{i,n}(\Theta _{n})}\) are constrained to eliminate arbitrage opportunities within this bond market, as dictated by the well-known Riccati equations.¹ For notational simplicity, we collect \(J\) bond yields into the vector \(\boldsymbol{Y_{i,t}}=[y_{i,t}^{(1)}, y_{i,t}^{(2)},...,y_{i,t}^{(J)}]^\top\), the \(J\) intercepts into \(\boldsymbol{A_X(\Theta_i)}=[a_{i,1}(\Theta _{1}), a_{i,2}(\Theta _{2}) ,...,a_{i,J}(\Theta _{J})]^\top\) \(\in \mathbb{R}^J\), and the \(N\) slope coefficients into a \(J \times N\) matrix \(B_X(\Theta_i)=[\boldsymbol{b_{i,1}(\Theta _{1})}^\top, \boldsymbol{b_{i,2}(\Theta _{2})}^\top, ...,\boldsymbol{b_{i,J}(\Theta _{J})}^\top]^\top\). Accordingly, the yield curve cross-section dimension of country \(i\) is: \[\begin{equation} \underset{(J \times 1)}{\vphantom{\Big|} \boldsymbol{Y_{i,t}}} = \underset{(J \times 1)}{\vphantom{\Big|} \boldsymbol{A_X(\Theta_i)}} + \underset{(J \times N)}{\vphantom{\Big|} B_X(\Theta_i)} \underset{(N \times 1)}{\vphantom{\Big|} \boldsymbol{X_{i,t}}}\text{.} \tag{4} \end{equation}\]

It follows from Equations (1) and (2) that the parameter set \(\Theta_i =\{\boldsymbol{\mu^Q_{i,X}},\Phi^Q_{i,X}, \Gamma_{i,X}, \delta_{i,0}, \boldsymbol{\delta_{i,1}}\}\) fully characterizes the cross-section of country’s \(i\) term structure. Importantly, Dai and Singleton (2000) demonstrate that this system is not identified without additional restrictions, since \(\boldsymbol{X_{i,t}}\) and any invertible affine transformation of \(\boldsymbol{X_{i,t}}\) yield observationally equivalent representations. To circumvent this problem, JPS (2014) adopt the three sets of (minimal) restrictions proposed by Joslin et al. (2011). First, they impose the latent factors to be zero-mean processes, forcing \(\boldsymbol{\mu^{Q}_{i,X}}= \boldsymbol{0}_N\). Second, they choose \(\boldsymbol{\delta_{i,1}}\) to be a \(N\)-dimensional vector whose entries are all equal to one. Lastly, \(\Phi^Q_{i,X}\) is a diagonal matrix, the elements of which are the real and distinct eigenvalues, \(\lambda^Q_i\), of the matrix of eigenvectors of \(\Phi^Q_{i,X}\). Based on this restriction set, Joslin et al. (2011) show that no additional invariant rotation is possible.

Joslin et al. (2011) also show that a rotation from \(\boldsymbol{X_{i,t}}\) to portfolios of yields, the spanned factors \(\boldsymbol{P_{i,t}}\), leads to an observationally equivalent model representation. This invariant transformation implies that \(N\) portfolios of yields are perfectly priced and observed without errors, while the remaining \(J-N\) portfolios are priced and observed imperfectly. Specifically, the spanned factors are computed as \(\boldsymbol{P_{i,t}}=V_i\boldsymbol{Y_{i,t}}\), for a full-rank matrix \(V_i\). Based on this definition, Equation (4) can be rearranged as an affine function of \(\boldsymbol{P_{i,t}}\)

\[\begin{equation} \underset{(J \times 1)}{\vphantom{\Big|} \boldsymbol{Y_{i,t}}}= \underset{(J \times 1)}{\vphantom{\Big|} \boldsymbol{A_P(\Theta_i)}}+ \underset{(J \times N)}{\vphantom{\Big|} B_P(\Theta_i)} \underset{(N \times 1)}{\vphantom{\Big|} \boldsymbol{P_{i,t}}}\text{.} \tag{5} \end{equation}\] where \(\boldsymbol{A_P(\Theta_i)}= \mathrm{I}_n - B_X(\Theta_i ) \left[ V_iB_X(\Theta_i \right] ^{-1}V_i \boldsymbol{A_X(\Theta_i)}\) and \(B_P(\Theta_i)=B_X(\Theta_i) \left[ V_iB_X(\Theta_i ) \right]^{-1}\).

The rotation from \(\boldsymbol{X_{i,t}}\) to \(\boldsymbol{P_{i,t}}\) is convenient for two key reasons. First, \(\boldsymbol{P_{i,t}}\) contains directly observable yield curve factors (unlike \(\boldsymbol{X_{i,t}}\)), with its \(N\) elements mapping to traditional yield curve components. For instance, for \(N=3\) and \(V_i\) being the weight matrix that results from a principal component analysis, the portfolios of yields \(\boldsymbol{P_{i,t}}\) are commonly referred to as the level, slope, and curvature factors (see Section 6.1). Second, it enables a convenient decomposition of the likelihood function, facilitating both estimation and the interpretation of model parameters.

Multicountry specifications (joint Q-dynamics model classes)

The cross-section multicountry extension is formed by stacking the country yields, spanned factors, and intercepts from Equation (5) into, respectively, \(\boldsymbol{Y_t}=[\boldsymbol{Y_{1,t}}^\top, \boldsymbol{Y_{2,t}}^\top, ...,\boldsymbol{Y_{C,t}}^\top]^\top\), \(\boldsymbol{P_t}=[\boldsymbol{P_{1,t}}^\top, \boldsymbol{P_{2,t}}^\top, ..., \boldsymbol{P_{C,t}}^\top]^\top\), and \(\boldsymbol{A_P(\Theta)}=[\boldsymbol{A_P^\top(\Theta_1)}, \boldsymbol{A_P^\top(\Theta_2)}, ..., \boldsymbol{A_P^\top(\Theta_C)}]^\top\), where \(C\) denotes the number of countries in this economic system. Additionally, we set \(B_{P}(\Theta)\) as block diagonal, \(B_P(\Theta)=B_P(\Theta_1) \oplus B_P(\Theta_2) \oplus \dots \oplus B_P(\Theta_C)\), where \(\oplus\) refers to the direct sum symbol. Accordingly,
\[\begin{equation} \underset{(CJ \times 1)}{\vphantom{\Big|} \boldsymbol{Y_{t}}} = \underset{(CJ \times 1)}{\vphantom{\Big|} \boldsymbol{A_{P}(\Theta)}} + \underset{(CJ \times CN)}{\vphantom{\Big|} B_{P}(\Theta)} \underset{(CN \times 1)}{\vphantom{\Big|} \boldsymbol{P_{t}}}\text{.} \tag{6} \end{equation}\]

Model time series dimension (P-dynamics)

In the modelling frameworks implemented in the MultiATSM package, the risk factor dynamics under the \(\mathbb{P}\)-measure must include at least \(N\) domestic spanned factors (\(\boldsymbol{P_{i,t}}\)) and \(M\) domestic unspanned factors (\(\boldsymbol{M_{i,t}}\)), and may optionally include \(G\) global unspanned factors (\(\boldsymbol{M_t^W}\)), depending on the specification. The dynamics of these risk factors evolves as either an unrestricted or a restricted VAR models. The unrestricted case corresponds to the JPS specification, while the restricted setup encompasses the GVAR and JLL frameworks.

It is worth stressing the role of unspanned factors in shaping the yield curve developments. Although these factors are absent in the cross-section dimension of the models, they influence the dynamics of the spanned factors which, in turn, affect directly bond yields.

JPS-based models

The country-specific state vector, \(\boldsymbol{Z_{i,t}}\), is formed from stacking the global and domestic (unspanned and spanned) risk factors: \(\boldsymbol{Z_{i,t}} = [\boldsymbol{M_t^{W^\top}}\), \(\boldsymbol{M_{i,t}}^\top\), \(\boldsymbol{P_{i,t}}^\top]^\top\). As such, \(\boldsymbol{Z_{i,t}}\) is a \(R\)-dimensional vector, where \(R =G + K\) and \(K = M + N\). In JPS-based setups, \(\boldsymbol{Z_{i,t}}\) follows a standard unrestricted Gaussian VAR(1): \[\begin{align} & \underset{(R \times 1)}{\vphantom{\Big|} \boldsymbol{Z_{i,t}}} = \underset{(R \times 1)}{\vphantom{\Big|} \boldsymbol{C_i^{\mathbb{P}}}} + \underset{(R \times R)}{\vphantom{\Big|} \Phi_i^{\mathbb{P}}} \underset{ (R \times 1)}{\vphantom{\Big|} \boldsymbol{Z_{i,t-1}}} + \underset{(R \times R)}{\vphantom{\Big|} \Gamma_i} \underset{(R \times 1)}{\vphantom{\Big|} \boldsymbol{\varepsilon_{Z,t}^{\mathbb{P}}}}\text{,} & \boldsymbol{\varepsilon_{Z,t}^{\mathbb{P}}} \sim {\mathcal{N}_R}(\boldsymbol{0}_R,\mathrm{I}_R)\text{,}% \tag{7} \end{align}\] where \(\boldsymbol{C_i^{\mathbb{P}}}\) denotes the vector of intercepts; \(\Phi_i^{\mathbb{P}}\), the feedback matrix; and \(\Gamma_i\), the Cholesky factor (a lower triangular matrix).

GVAR-based models

In the MultiATSM package, the GVAR setup is formed from two parts: the marginal and the VARX\(^{*}\) models. The former captures the joint dynamics of the global economy, whereas the latter describes the developments from the domestic factors. For a thorough description of GVAR models, See Chudik and Pesaran (2016).

The marginal model is an unrestricted VAR(\(1\)) featuring exclusively the global factors: \[\begin{align} & \underset{(G \times 1)}{\vphantom{\Big|} \boldsymbol{M_t^W}}= \underset{(G \times 1)}{\vphantom{\Big|} \boldsymbol{C^W}} + \underset{(G \times G)}{\vphantom{\Big|} \Phi^W} \underset{(G \times 1)}{\vphantom{\Big|} \boldsymbol{M_{t-1}^W}} + \underset{(G \times G)}{\vphantom{\Big|} \Gamma^W}\underset{(G \times 1)}{\vphantom{\Big|} \boldsymbol{\varepsilon_{t}^W}}\text{,} & \boldsymbol{\varepsilon_t^W} \sim {\mathcal{N}_G}(\boldsymbol{0}_G,\mathrm{I}_G). \tag{8} \end{align}\]

The VARX\(^{*}\) setups are country-specific small-scale VAR models containing global factors and weakly exogenous ‘star’ variables — weighted averages of foreign variables — built as \[\begin{equation} \boldsymbol{Z_{i,t}^{\ast^\top}} = \sum_{j=1}^{C} w_{i,j} \boldsymbol{Z_{j,t}^\top}, \qquad \sum_{j=1}^{C} w_{i,j}= 1, \quad w_{i,i}=0 \quad \forall i \in \{1,2, ...C \}, \tag{9} \end{equation}\] where \(Z_{j,t}\) is a \(K-\)dimension vector of domestic factors \(\boldsymbol{Z_{j,t}} = [\boldsymbol{M_{j,t}}^\top\), \(\boldsymbol{P_{j,t}}^\top]^\top\) and \(w_{i,j}\) is a scalar that measures the degree of connectedness of country \(i\) with country \(j\).

These models follow a VARX\(^{*}(p,q,r)\) specification, where \(p\), \(q\) and \(r\) are the number of lags from, respectively, the domestic, the star, and the global risk factors. The MultiATSM package provides the estimates for the case \(p=q=r=1\). In such a case, the dynamics of \(\boldsymbol{Z_{i,t}}\) is described as a VARX\(^{*}\) of the following form: \[\begin{align} & \underset{(K \times 1)}{\vphantom{\Big|} \boldsymbol{Z_{i,t}}} = \underset{(K \times 1)}{\vphantom{\Big|} \boldsymbol{C^X_{i}}} + \underset{(K \times K)}{\vphantom{\Big|} \Phi^X_{i}} \underset{(K \times 1)}{\vphantom{\Big|} \boldsymbol{Z_{i,t-1}}} + \underset{(K \times K)}{\vphantom{\Big|} \Phi^{X^\ast}_i} \underset{(K \times 1)}{\vphantom{\Big|} \boldsymbol{Z_{i,t-1}^{\ast}}} + \underset{(K \times G)}{\vphantom{\Big|} \Phi_{i}^{X^{W}}} \underset{(G \times 1)}{\vphantom{\Big|} \boldsymbol{M_{t-1}^{W}}} + \underset{(K \times K)}{\vphantom{\Big|} \Gamma_{i}^{X}} \underset{(K \times 1)}{\vphantom{\Big|} \boldsymbol{\varepsilon^X_{i,t}}}\text{,} & \boldsymbol{\varepsilon^X_{i,t}} \sim {\mathcal{N}_K}(\boldsymbol{0}_K,\mathrm{I}_K). \tag{10} \end{align}\]

Additionally, GVAR models require, as an intermediate step, the specification of country-specific \(2K \times CK\)-link matrices, \(W_i\), to unify the individual VARX\(^{*}\) models. Formally, \[\begin{equation} \begin{bmatrix} \boldsymbol{Z_{i,t}} \\ \boldsymbol{Z_{i,t}}^{*} \end{bmatrix}_{2K \times 1} \equiv \underset{(2K \times CK)}{W_i} \begin{bmatrix} \boldsymbol{Z_{1,t}} \\ \boldsymbol{Z_{2,t}} \\ \vdots \\ \boldsymbol{Z_{C,t}} \end{bmatrix}_{CK \times 1}. \tag{11} \end{equation}\]

Last, to compose the \(F\)-dimensional state vector for \(F = G + CK\), we gather the global economic variables and the country-specific risk factors, as \(\boldsymbol{Z_t} = [\boldsymbol{M_{t}^{W^\top}}\), \(\boldsymbol{Z_{1,t}}^\top\), \(\boldsymbol{Z_{2,t}}^\top, \hdots \boldsymbol{Z_{C,t}}^\top]^\top\). As such, we can form a first order GVAR process as \[\begin{align} & \underset{(F \times 1)}{\vphantom{\Big|} \boldsymbol{Z_t}} = \underset{(F \times 1)}{\vphantom{\Big|} \boldsymbol{C_y}} + \underset{(F \times F)}{\vphantom{\Big|} \Phi_y} \underset{(F \times F)}{\vphantom{\Big|} \boldsymbol{Z_{t-1}}} + \underset{(F \times F)}{\vphantom{\Big|} \Gamma_y} \underset{(F \times 1)}{\vphantom{\Big|} \boldsymbol{\varepsilon_{y,t}}}\text{,} & \boldsymbol{\varepsilon_{y,t}} \sim {\mathcal{N}_F}(\boldsymbol{0}_F,\mathrm{I}_F)\text{,} \tag{12} \end{align}\] where \(\boldsymbol{C_y} = [\boldsymbol{C^{W^\top}}\), \(\boldsymbol{C_1^{X^\top}}\), \(\boldsymbol{C_2^{X^\top}}\),… \(\boldsymbol{C_C^{X^\top}}]^\top\), \(\boldsymbol{\varepsilon_{y,t}} =[ \boldsymbol{\varepsilon^{W^\top}_t}\), \(\boldsymbol{\varepsilon_{1,t}^{X^\top}}\), \(\boldsymbol{\varepsilon_{2,t}^{X^\top}}\) … \(\boldsymbol{\varepsilon_{C,t}^{X^\top}}]^\top\), \(\Gamma_y=\Gamma^W \oplus \Gamma_1^X \oplus \Gamma_2^X \oplus \dots \oplus \Gamma_C^X\), and

\[\begin{equation} \Phi_y = \begin{bmatrix} \Phi^W & 0_{\scriptscriptstyle{G \times CK}} \\ \Phi^{X^{W}} & G_1 \end{bmatrix}_{F \times F} , \end{equation}\]

where \(\Phi^{X^{W}}= \begin{bmatrix} \Phi^{X^{W}}_1 \\ \Phi^{X^{W}}_2 \\ \vdots \\ \Phi^{X^{W}}_C \end{bmatrix}_{CK \times G}\) and \(G_1= \begin{bmatrix} \Phi_1W_1 \\ \Phi_2W_2 \\ \vdots \\ \Phi_CW_C \end{bmatrix}_{CK \times CK}\), for \(\Phi_i= [\Phi_i^{X}\), \(\Phi_i^{X^*}]\) and \(\quad \forall i \in \{1,2, ...C \}\).

JLL-based models

JLL-based models incorporate three components: (i) the global economy, (ii) a dominant large economy,² and (iii) a set of smaller economies. The state vector is formed from a number of linear projections to build domestic risk factors that are free of the influence of the variables from other countries and/or from the global economy.

The construction of the domestic spanned factors proceeds in two steps. First, for each economy \(i\), \(\boldsymbol{P_{i,t}}\) is projected on \(\boldsymbol{M_{i,t}}\) of this same country \[\begin{equation} \underset{ (N \times 1)}{\vphantom{\Big|} \boldsymbol{P_{i,t}}} = \underset{(N \times M)}{\vphantom{\Big|} b_i} \underset{(M \times 1)}{\vphantom{\Big|} \boldsymbol{M_{i,t}}} + \underset{ (N \times 1)}{\vphantom{\Big|} \boldsymbol{P_{i,t}^e}} \text{,} \tag{13} \end{equation}\] where the residuals \(\boldsymbol{P_{i,t}^e}\) are orthogonal to the economic fundamentals of the country \(i\).

Second, for the non-dominant economies, \(\boldsymbol{P_{i,t}^e}\) is additionally projected on the orthogonalized spanned factors of the dominant country, indexed by \(D\), as follows: \[\begin{equation} \underset{(N \times 1)}{\vphantom{\Big|} \boldsymbol{P_{i,t}^e}} = \underset{(N \times N)}{\vphantom{\Big|} c_i^D} \underset{(N \times 1)}{\vphantom{\Big|} \boldsymbol{P_{D,t}^e}} + \underset{(N \times 1)}{\vphantom{\Big|} \boldsymbol{P_{i,t}^{e*}}}\text{,} \quad i \neq D, \tag{14} \end{equation}\] where \(\boldsymbol{P_{i,t}^{e*}}\) corresponds to the non-dominant country \(i\) set of residuals.

The design of the domestic unspanned factors also features two steps: for the dominant economy, \(\boldsymbol{M_{D,t}}\) is projected on the global economic factors \[\begin{equation} \underset{(M \times 1)}{\vphantom{\Big|} \boldsymbol{M_{D,t}}} = \underset{(M \times G)}{\vphantom{\Big|} a_D^W} \underset{(G \times 1)}{\vphantom{\Big|} \boldsymbol{M_t^W}} + \underset{(M \times 1)}{\vphantom{\Big|} \boldsymbol{M_{D,t}^e}} \text{,} \tag{15} \end{equation}\] and, for the other economies, the residuals of the previous regression are used to compute \[\begin{equation} \underset{(M \times 1)}{\vphantom{\Big|} \boldsymbol{M_{i,t}}} = \underset{(M \times G)}{\vphantom{\Big|} a_i^W} \underset{(G \times 1)}{\vphantom{\Big|} \boldsymbol{M_t^W}} + \underset{(M \times M)}{\vphantom{\Big|} a_i^D} \underset{(M \times 1)}{\vphantom{\Big|} \boldsymbol{M_{D,t}^e}} + \underset{(M \times 1)}{\vphantom{\Big|} \boldsymbol{M_{i,t}^{e*}}}\text{.} \tag{16} \end{equation}\]

Accordingly, the state vector is formed by \(\boldsymbol{Z_t^e}= [\boldsymbol{M_t^{W^\top}}\), \(\boldsymbol{M_{D,t}^{e^\top}}\), \(\boldsymbol{P_{D,t}^{e^\top}}\), \(\boldsymbol{M_{2,t}^{e*^\top}}\), \(\boldsymbol{P_{2,t}^{e*^\top}}\) … \(\boldsymbol{M_{C,t}^{e*^\top}}\), \(\boldsymbol{P_{C,t}^{e*^\top}}]^\top\) and its dynamics evolve as a restricted VAR(1), \[\begin{align} & \underset{(F \times 1)}{\vphantom{\Big|} \boldsymbol{Z_t^e}}= \underset{(F \times 1)}{\vphantom{\Big|} \boldsymbol{C^{e}_Y}} + \underset{(F \times F)}{\vphantom{\Big|} \Phi^e_Y} \underset{(F \times 1)}{\vphantom{\Big|} \boldsymbol{Z_{t-1}^e}} + \underset{(F \times F)}{\vphantom{\Big|} \Gamma_{Y}^e} \underset{(F \times 1)}{\vphantom{\Big|} \boldsymbol{\varepsilon^e_{Z,t}}}\text{,} & \boldsymbol{\varepsilon _{Z,t}^{e}} \sim {\mathcal{N}_F}(\boldsymbol{0}_F,\mathrm{I}_F). \tag{17} \end{align}\] JLL (2015) impose a set of zero restrictions on \(\Phi^e_Y\) and \(\Gamma_{Y}^e\), with their detailed structure provided in the original study.

Estimation procedures

The approach proposed by JPS (2014) enables an efficient estimation procedure through its structural design. Specifically, the parameters governing the \(\mathbb{Q}\)- and \(\mathbb{P}\)-measures can be estimated independently. The only exception is the variance-covariance matrix, \(\Sigma\), which appears in both likelihood functions and, therefore, must be estimated jointly.

In JLL (2015), however, the authors adopt a simplified estimation procedure by estimating the \(\Sigma\) matrix exclusively under the \(\mathbb{P}\)-measure. While they acknowledge that this approach is not fully efficient, they argue that the empirical implications are limited in their application.

3 The ATSMs available at the MultiATSM package

As outlined in the previous section, the ATSMs implemented in the MultiATSM package differ in the specification of their \(\mathbb{Q}\)- and \(\mathbb{P}\)-measure dynamics. In short, under the \(\mathbb{Q}\)-measure, models can be specified either on a country-by-country basis (JPS, 2014) or jointly across countries (JLL, 2015; CM, 2024). Under the \(\mathbb{P}\)-measure, risk factor dynamics follow a VAR(1) process, which may be unrestricted, as in the JPS-related frameworks, or restricted, as in the JLL and GVAR specifications.

MultiATSM provides support for eight different classes of ATSMs based on these modelling approaches. These classes vary along several dimensions: the specification of the \(\mathbb{P}\)- and \(\mathbb{Q}\)-dynamics, the estimation approach, and whether a dominant economy is included. Table 1 summarizes the defining features of each model class available in the package. A brief overview of these specifications follows below.

Table 1: Summary of model features

		P-dynamics							Q-dynamics				Sigma matrix estimation				Dom. Eco.
		Single		Joint					Single	Joint			P only	P and Q
		UR	R	UR	R
					JLL	GVAR
Unrestricted VAR
JPS original		x							x					x
JPS global				x					x					x
JPS multi				x						x				x
Restricted VAR (GVAR)
GVAR single						x			x					x
GVAR multi						x				x				x
Restricted VAR (JLL)
JLL original					x					x			x				x
JLL No DomUnit					x					x			x
JLL joint Sigma					x					x				x			x
Note:
Risk factor dynamics under the \(\mathbb{P}\)-measure may follow either an unrestricted (UR) or a restricted (R) specification. The set of restrictions present in the JLL-based and GVAR-based models are described in Jotikasthira et al. (2015) and Candelon and Moura (2024), respectively. The estimation of the \(\Sigma\) matrix is done either exclusively with the other parameters of the \(\mathbb{P}\)-dynamics (P column) or jointly under both \(\mathbb{P}\)- and \(\mathbb{Q}\)-parameters (P and Q column). Dom. Eco. relates to the presence of a dominant economy. The entries featuring x indicate that the referred characteristic is part of the model.

The ATSMs in which the estimation is performed separately for each country are labeled as JPS original, JPS global and GVAR single. In the JPS original setup, the set of risk factors includes exclusively each country’s domestic variables and the global unspanned factors, whereas JPS global and GVAR single also incorporate domestic risk factors of the other countries of the economic system. Noticeably, the difference between JPS global and GVAR single stem from the set of restrictions imposed under the \(\mathbb{P}\)-dynamics.

Within the multicountry frameworks, certain features are worth noting. The JLL original model reproduces the setup in JLL (2015), assuming an economic cohort composed of a globally dominant economy and a set of smaller countries, and estimating the \(\Sigma\) matrix exclusively under the \(\mathbb{P}\)-measure. The two alternative versions assume the absence of a dominant country (JLL No DomUnit) and the estimation of \(\Sigma\) under both the \(\mathbb{P}\) and \(\mathbb{Q}\) measures (JLL joint Sigma), as in JPS (2014). The remaining specifications differ in their \(\mathbb{P}\)-dynamics: either by an unrestricted VAR model (JPS multi) or by a GVAR setup (GVAR multi), as proposed in CM (2024).

4 Package dataset

The MultiATSM package provides datasets that approximate those used in the GVAR-based ATSMs of Candelon and Moura (2023) and CM (2024). The data requirements for estimating GVAR models encompass those of all other model classes, making them suitable for generating outputs across all models supported by the package. As such, the examples in the following sections use the dataset from CM (2024).

The LoadData() function provides access to the datasets included in the package. To load the data from CM (2024), set the argument to CM_2024:

LoadData("CM_2024")

This function returns three sets of data. The first contains time series of zero-coupon bond yields for four emerging market economies: China, Brazil, Mexico, and Uruguay. The data spans monthly intervals from June 2004 to January 2020. For the purpose of model estimation, the package requires that (i) bond yield maturities are the same across all countries;³ and (ii) yields must be expressed in annualized percentage terms (not basis points). Note that the MultiATSM package does not provide routines for bootstrapping zero-coupon yields from coupon bonds, so any such treatment must be handled by the user.

The second dataset comprises time series for unspanned risk factors — specifically, the macroeconomic indicators economic growth and inflation — covering the same period as the bond yield data. These data cover both (i) domestic variables for each of the four countries in the sample and (ii) corresponding global indicators. The construction of unspanned risk factors, like that of bond yields, must be carried out externally by the user.

The final dataset contains measures of interconnectedness, proxied by trade flows, which are specifically required for estimating the GVAR-based models. The trade flow data report the annual value of goods imported and exported between each pair of countries in the sample, starting from 1948. All values are expressed in U.S. dollars on a free-on-board basis. These data are used to construct the transition matrix in the GVAR framework.

5 Required user inputs

Fundamental inputs

To estimate any model, the user must specify several general inputs, which can be grouped into the following categories:

1. Desired ATSM class (ModelType): a character vector containing the label of the model to be estimated as described in Table 1;

2. Risk Factor Features. This includes the following list of elements:

• Set of economies (Economies): a character vector containing the names of the economies which are part of the economic system;

• Global variables (GlobalVar): a character vector containing the labels of the \(G\) global unspanned factors. Studies examining the impact of global developments on bond prices could include proxy measures of global inflation and global economic activity in this category (Jotikasthira et al. 2015; Abbritti et al. 2018; Candelon and Moura 2024);

• Domestic variables (DomVar): a character vector containing the labels of the \(M\) domestic unspanned factors. These typically correspond to measures of domestic inflation and economic activity, the standard macroeconomic indicators monitored by central banks (Ang and Piazzesi 2003; Joslin et al. 2014; Jotikasthira et al. 2015; Candelon and Moura 2024);

• Number of spanned factors (\(N\)): a scalar representing the number of country-specific spanned factors. Although, in principle, \(N\) could vary across countries, the models provided in the package assume a common value of \(N\) for all countries. A common choice in the literature is \(N=3\), as in JPS (2014) and CM (2024), since this produces an excellent cross-sectional fit of bond yields (Litterman and Scheinkman (1991)). Other studies, such as Adrian et al. (2013), extend the specification to \(N=5\), arguing that it improves the performance of model-implied term premia. Further intuition on the role and interpretation of spanned factors is provided in Section 6.1.

3. Sample span:

• Initial sample date (t0): the start of the sample period in the format dd-mm-yyyy;

• End sample date (tF): the end of the sample period in the format dd-mm-yyyy.

4. Data Frequency (DataFreq): a character vector specifying the frequency of the time series data. The available options are: Annually, Quarterly, Monthly, Weekly, Daily Business Days, and Daily All Days;

5. Stationarity constraint under the \(\mathbb{Q}\)-dynamics (StatQ): a logical that takes TRUE if the user wishes to impose that the largest eigenvalue under the \(\mathbb{Q}\)-measure, \(\lambda^Q_i\), is strictly less than 1. While enforcing this stationarity constraint may increase estimation time, it can improve convergence and numerical stability. Moreover, by inducing near-cointegration, the eigenvalue restriction helps to pin down more plausible dynamics for bond risk premia (Bauer et al. 2012; Joslin et al. 2014);

6. Selected folder to save the graphical outputs (Folder2Save): path where the selected graphical outputs will be saved. If set to NULL, the outputs are stored in the user’s temporary directory (accessible via tempdir());

7. Output label (OutputLabel): A single-element character vector containing the name used in the file name that stores the model outputs.

The following provides an example of the basic model input specification:

ModelType <- "JPS original"
Economies <- c("Brazil", "Mexico", "Uruguay")
GlobalVar <- c("Gl_Eco_Act", "Gl_Inflation")
DomVar <- c("Eco_Act", "Inflation")
N <- 3
t0 <- "01-07-2005"
tF <- "01-12-2019"
DataFreq <- "Monthly"
StatQ <- FALSE
Folder2Save <- NULL
OutputLabel <- "Model_demo"

Model-specific inputs

GVARlist and JLLlist

The inputs described above are sufficient for estimating all variants of the JPS models presented in Table 1. However, estimating the GVAR or JLL setups requires additional elements. For clarity, these extra inputs should be organized into separate lists for each model. This section outlines the general structure of both lists, while Section 6.2 provide a more detailed explanations of their components and available options, reflecting the broader scope of each setup.

For GVAR models, the required inputs are twofold. First, the user must specify the dynamic structure of each country’s VARX model. For example:

VARXtype <- "unconstrained"

Next, provide the desired inputs to build the transition matrix. For instance:

data('TradeFlows')
W_type <- "Sample Mean"
t_First_Wgvar <- "2000"
t_Last_Wgvar <- "2015"
DataConnectedness <- TradeFlows 

Based on these inputs, a complete instance of the GVARlist object is

GVARlist <- list(VARXtype = "unconstrained", W_type = "Sample Mean", 
                 t_First_Wgvar = "2000", t_Last_Wgvar = "2015", 
                 DataConnectedness = TradeFlows) 

For the JLL frameworks, if the chosen model is either JLL original or JLL joint Sigma, it suffices to specify the name of the dominant economy. Otherwise, for the JLL No DomUnit class, the user must set None. For instance:

## Example for "JLL original" and "JLL joint Sigma" models
JLLlist <- list(DomUnit =  "China")

## For "JLL No DomUnit" model
JLLlist <- list(DomUnit =  "None")

BRWlist

In an influential paper, Bauer et al. (2012) (henceforth BRW, 2012) show that estimates from traditional ATSMs often suffer from severe small-sample bias. This can lead to unrealistically stable expectations for future short-term interest rates and, consequently, distort term premium estimates for long-maturity bonds. To address this issue, BRW (2012) propose an indirect inference estimator based on a stochastic approximation algorithm, which corrects for bias and enhances the persistence of short-term interest rates, resulting in more plausible term premium dynamics.

It is worth noting that this framework serves as a complementary feature to the core ATSMs and can therefore be applied to any of the model types supported by the MultiATSM package. If the user intends to implement a model following the BRW (2012) approach, a few additional inputs must be specified. These include:

• Mean or median of physical dynamic estimates (Cent_Measure): compute the mean or the median of the \(\mathbb{P}\)-dynamics estimates after each bootstrap iteration by setting the option to Mean (for the mean) or Median (for the median);

• Adjustment parameter (gamma): this parameter controls the degree of shrinkage applied to the difference between the estimates prior to the bias correction and the bootstrap-based estimates after each iteration. It remains fixed across iterations and must lie in the interval \((0,1)\);

• Number of iteration (N_iter) : total number of iterations used in the stochastic approximation algorithm after burn-in;

• Number of bootstrap samples (B): quantity of simulated samples used in each burn-in or actual iteration;

• Perform closeness check (checkBRW): indicates whether the user wishes to compute the root mean square distance between the model estimates obtained through the bias-correction method and those generated via no bias-correction. The default is TRUE;

• Number of bootstrap samples used in the closeness check (B_check): default is equal to 100,000 samples;

• Eigenvalue restriction (Eigen_rest): impose a restriction on the largest eigenvalue under the \(P\)-measure after applying the bias correction procedure. Default is \(1\);

• Number of burn-in iteration (N_burn): quantity of the iterations to be discarded in the first stage of the bias-correction estimation process. The recommended number is \(15\%\) of the total number of iterations. In practice, this resembles the burn-in concept in Markov chain Monte Carlo methods. Particularly, the BRW (2012) stochastic approximation algorithm is iterative, and for a sufficiently large number of iterations, the parameters converge to their true values. As such, discarding early iterations avoids the need for assessing a computationally costly exit condition.

BRWlist <- within(list(Cent_Measure = "Mean", gamma = 0.2, N_iter = 500, B = 50, 
                       checkBRW = TRUE, B_check = 1000, Eigen_rest = 1), 
                       N_burn <- round(N_iter * 0.15))

Additional inputs for numerical and graphical outputs

Once the desired features are selected and the parameters of the chosen ATSM have been estimated, the MultiATSM package provides tools to generate the following numerical and graphical outputs via the NumOutputs() function:

• Time-series dynamics of the risk factors;
• Model fit of the bond yields;
• Orthogonalized impulse response functions (IRFs);
• Orthogonalized forecast error variance decompositions (FEVDs);
• Generalized impulse response functions (GIRFs);
• Generalized forecast error variance decompositions (GFEVDs);
• Decomposition of bond yields into expected and term premia components.

These outputs are organized into distinct analytical components, each offering different insights into model behavior and its economic interpretation.

The time-series dynamics of the risk factors are displayed in separate subplots: one for each global factor, and one subplot per domestic risk factor showing all countries in the economic system. The model fit of the bond yields is provided through two measures of model-implied yields. The first is a fitted measure derived solely from the cross-sectional component, as in Equation (5) for single-country models and Equation (6) for multicountry setups. This measure reflects the fit based exclusively on the parameters governing the \(\mathbb{Q}\)-dynamics. The second incorporates both the physical and risk-neutral dynamics, combining the cross-sectional equations with the state evolution specified by each ATSM.

The impulse response functions and variance decompositions are available in both orthogonalized and generalized forms. The orthogonalized outputs (IRFs and FEVDs) are computed using a short-run recursive identification scheme, meaning they depend on the ordering of the selected risk factors. Specifically, the package is structured to place global unspanned factors first, followed by its domestic unspanned and spanned factors within each country, in the order in which countries are listed in the Economies vector. In contrast, the generalized versions (GIRFs and GFEVDs) are robust to factor ordering but allow for correlated shocks across risk factors (Pesaran and Shin 1998). For the numerical computation of these outputs, a horizon of analysis has to be specified, e.g., Horiz <- 100.

The bond yield decomposition can be performed with respect to two measures of risk compensation: term premia and forward premia. While the term premium is derived directly from the bond yield levels, the forward premium is obtained from the decomposition of forward rates. A more formal presentation of both measures is provided in the Appendix.

Users must specify the desired graph types in a character vector. Available options include: RiskFactors, Fit, IRF, FEVD, GIRF, GFEVD, and TermPremia. For example:

DesiredGraphs <- c("Fit", "GIRF", "GFEVD", "TermPremia")

Moreover, for all models, users must indicate the types of variables of interest (yields, risk factors, or both). For JLL-type models specifically, users must also specify whether to include the orthogonalized versions. Each of these options should be set to TRUE to generate the corresponding graphs, and FALSE otherwise.

WishGraphRiskFac <- FALSE
WishGraphYields <- TRUE
WishOrthoJLLgraphs <- FALSE

The desired graphical outputs are stored in the selected folder, Folder2Save. Alternatively, users can display the desired plots directly in the console without saving them to Folder2Save by using the autoplot() method.

Bootstrap settings

Horowitz (2019) shows that bootstrap methods generally produce more accurate statistical inference than those based on asymptotic distribution theory. To generate confidence intervals using bootstrap, via the Bootstrap() function, an additional list of inputs must be provided:

• Desired bootstrap procedure (methodBS): the user must select one of the following options: (i) standard residual bootstrap (bs); (ii) wild bootstrap (wild); or (iii) block bootstrap (block). If the block bootstrap is selected, the block length must also be specified. The residual bootstrap is a conventional method that is straightforward to implement when a parametric model, such as a VAR model, is available. The block bootstrap makes weaker assumptions about the data-generating process and is well-suited to handling both weak and strong serial dependence. The wild bootstrap is particularly appropriate for data exhibiting heteroskedasticity (Horowitz 2019);
• Number of bootstrap draws (ndraws): Kilian and Lütkepohl (2017) suggest that, in VAR specifications, ndraws can range from a few hundred to several thousand, depending on factors such as sample size, lag order, and the desired quantiles of the distribution. Illustrating this, CM (2024) set \(ndraws= 1,000\) in their ATSM to construct confidence intervals for IRFs;
• Confidence level expressed (pctg): the desired confidence level should be expressed in percentage points. Common choices in VAR-related setups include \(68\%\), \(90\%\) and \(95\%\) (Kilian and Lütkepohl 2017).

Bootlist <- list(methodBS = 'block', BlockLength = 4, ndraws =  1000, pctg   =  95)

Out-of-sample forecast settings

To generate bond yield forecasts, use ForecastYields() with the following inputs:

• Forecast horizon (ForHoriz): Number of forecast horizons in periods;
• Index of the first observation (t0Sample): time index of the first observation included in the information set;
• Index of the last observation (t0Forecast): time index of the last observation in the information set used to generate the first forecast;
• Method used for forecast computation (ForType): forecasts can be generated using either a rolling or expanding window. To use a rolling window, set this parameter to Rolling. In this case, the sample length for each forecast is fixed and defined by t0Sample. For expanding window forecasts, set this input to Expanding, allowing the information set to increase at each forecast iteration.

ForecastList <- list(ForHoriz = 12, t0Sample = 1, t0Forecast = 70, ForType = "Rolling")

6 Model estimation

Using the dataset described in Section 4, the estimation of the ATSM proceeds in three main steps. First, the country-specific spanned factors are estimated, which, along with the global and domestic unspanned factors, form the complete set of risk factors used in the subsequent estimation steps. Second, the package estimates the parameters governing the dynamics of the risk factors under the \(\mathbb{P}\)-measure. Finally, it optimizes the full ATSM specification, including the parameters under the \(\mathbb{Q}\)-measure.

As will be made clear in Section 7, although the functions introduced in this section can be used individually, they are primarily designed to be used together with the broader set of functions available in the MultiATSM package. However, as these functions play a central role in the package structure, they warrant a dedicated section.

Spanned factors

The spanned factors for country \(i\), denoted by \(\boldsymbol{P_{i,t}}\), are typically obtained as the first \(N\) principal components (PCs) of the observed bond yields. The PC method provides orthogonal linear combinations of the original variables, ordered by their ability to capture the variance in the data. Formally, \(\boldsymbol{P_{i,t}}\) is computed as \(\boldsymbol{P_{i,t}} = w_i \boldsymbol{Y_{i,t}}\), where yields are ordered by increasing maturity in \(\boldsymbol{Y_{i,t}}\), and \(w_i\) is the matrix of eigenvectors derived from the covariance matrix of \(\boldsymbol{Y_{i,t}}\).

In the case of \(N = 3\), the spanned factors are traditionally interpreted as level, slope, and curvature components of the yield curve (Litterman and Scheinkman 1991). This interpretation stems from the properties of the \(w_i\) matrix, as illustrated below:

data('Yields')
w <- pca_weights_one_country(Yields, Economy = "Uruguay")

In matrix w, each row holds the weights for constructing a spanned factor. The first row relates to the level factor, with weights loading roughly equally across maturities. As such, high (low) values of the level factor indicate an overall high (low) value of yields across all maturities. The second row features increasing weights with maturity, capturing the slope of the yield curve: high values indicate steep curves, while low values reflect flat or inverted curves. The third row corresponds to the curvature factor, with weights emphasizing medium-term maturities. This captures the ‘hump-shaped’ features of the yield curve typically associated with changes in its curvature. These concepts are also graphically illustrated in Figure 1.

Figure 1: Yield loadings on the spanned factors. Example using bond yield data for Uruguay. Graph generated using the ggplot2 package (Wickham 2016).

The user can directly obtain the time series of the country-specific spanned factors by calling Spanned_Factors(), as shown below:

data('Yields')
Economies <- c("China", "Brazil", "Mexico", "Uruguay")
N <- 2
SpaFact <- Spanned_Factors(Yields, Economies, N)

The P-dynamics estimation

As presented in Table 1 and explained in detail in Section 2, the dynamics of the risk factors under the \(\mathbb{P}\)-measure in the available models follow a VAR(1) process. This specification can be fully unrestricted, as in the JPS-related models, or subject to restrictions, as in the GVAR and JLL frameworks. This subsection illustrates how each of these model configurations is implemented.

VAR

To use VAR(), the user needs to select the appropriate set of risk factors for the model being estimated and specify unconstrained in the argument VARtype. In the two examples presented below, the outputs are the intercept vector, the feedback matrix, and the variance–covariance matrix for a VAR(1) model under the \(\mathbb{P}\)-measure:

## Example 1: "JPS global" and "JPS multi" models
data("RiskFacFull")
PdynPara <- VAR(RiskFacFull, VARtype = "unconstrained")

## Example 2: "JPS original" model for China 
FactorsChina <- RiskFacFull[1:7, ]
PdynPara <- VAR(FactorsChina, VARtype = "unconstrained")

GVAR

The GVAR() function estimates a GVAR(1) model constructed from country-specific VARX\(^{*}(1,1,1)\) specifications. It requires two main inputs: the number of domestic spanned factors (\(N\)) and a set of elements grouped in the GVARinputs list. The latter consists of four components:

1. Economies: a \(C-\)dimensional character vector containing the names of the economies present in the economic system;

2. GVAR list of risk factors: a list of risk factors sorted by country in addition to the global variables. An example of the expected data structure is:

data("GVARFactors")

To assist in formatting the data accordingly, users may use the DatabasePrep() function;

3. VARX type: a character vector specifying the desired structure of the VARX\(^{*}\) model. Two general options are available:

• Fully unconstrained: specify as unconstrained. This option estimates each equation in the system separately via ordinary least squares, without imposing any restrictions.

• With constraints: imposes specific set of zero restrictions on the feedback matrix. This category includes two sub-options: (a) constrained: Spanned Factors prevents foreign spanned factors from affecting any domestic risk factor; (b) constrained: [factor name] restricts the specified risk factor to be influenced only by its own lags and the lags of its associated star variables. In both cases, the VARX\(^{*}\) is estimated using restricted least squares.

data('GVARFactors')
GVARinputs <- list(Economies = Economies, GVARFactors = GVARFactors, 
                   VARXtype ="constrained: Inflation")

4. Transition matrix: a \(C \times C\) matrix that captures the degree of interdependence across the countries in the system. Each entry \((i,j)\) represents the strength of the dependence of economy \(i\) on economy \(j\). As an example, the matrix below is computed from bilateral trade flow data, averaged over the period 2006–2019, for a system comprising China, Brazil, Mexico, and Uruguay. The rows are normalized so that the weights sum to \(1\) for each country (i.e., each row of the matrix sums to \(1\)). The transition matrix can be generated using Transition_Matrix(), as illustrated in the Appendix :

         China Brazil Mexico Uruguay
China   0.0000 0.6549 0.3155  0.0296
Brazil  0.8269 0.0000 0.1234  0.0497
Mexico  0.8596 0.1326 0.0000  0.0078
Uruguay 0.3811 0.5498 0.0691  0.0000

With inputs specified, the user can estimate a GVAR model using:

data("GVARFactors")
GVARinputs <- list(Economies = Economies, GVARFactors = GVARFactors,  
                   VARXtype = "unconstrained", Wgvar = W_gvar)
N <- 3
GVARpara <- GVAR(GVARinputs, N, CheckInputs = TRUE)

Note that the CheckInputs parameter should be set to TRUE to perform a consistency check on the inputs specified in GVARinputs prior to the \(\mathbb{P}\)-dynamics estimation.

JLL

The JLL() function estimates the physical parameters. Required inputs are:

1. Risk factors: a time series matrix of the risk factors in their non-orthogonalized form;

2. Number of spanned factors (\(N\)): a scalar representing the number of country-specific spanned factors;

3. JLLinputs: a list object containing the following elements:

• Economies: a \(C\)-dimensional character vector listing the economies;

• Dominant Economy: a character vector indicating either the name of the country assigned as the dominant economy (for JLL original and JLL jointSigma models), or None (for the JLL No DomUnit case);

• Estimate Sigma Matrices: a logical equal to TRUE if the user wishes to estimate the full set of JLL sigma matrices (i.e., variance-covariance and Cholesky factor matrices), and FALSE otherwise. Since this numerical estimation is costly, it may significantly increase computation time;

• Precomputed Variance-Covariance Matrix: in some instances, a precomputed variance-covariance matrix from the non-orthogonalized dynamics can be supplied here to save time and memory. If no such matrix is available, this input should be set to NULL;

• JLL type: a character string specifying the chosen JLL model, following the classification described in Table 1.

## First set the JLLinputs 
ModelType <- "JLL original"   
JLLinputs <- list(Economies = Economies, DomUnit = "China", WishSigmas = TRUE,  
                  SigmaNonOrtho = NULL, JLLModelType = ModelType)

## Then, estimate the desired the P-dynamics from the desired JLL model
data("RiskFacFull")
N <- 3
JLLpara <- JLL(RiskFacFull, N, JLLinputs, CheckInputs = TRUE)

The CheckInputs input is set to TRUE to perform a consistency check on the inputs specified in JLLinputs before running the \(\mathbb{P}\)-dynamics estimation.

ATSM estimation

Estimating the ATSM parameters involves maximizing the log-likelihood function to obtain the best-fitting model parameters using Optimization(). The unspanned risk factor framework of JPS (2014) (and, therefore, all its multicountry extensions) follows a model parameterization similar to that proposed in Joslin et al. (2011). Particularly, it requires estimating a set of six parameter blocks:

1. The risk-neutral long-run mean of the short rate (\(r0\));

2. The risk-neutral feedback matrix (\(K1XQ\));

3. Standard deviation of measurement errors for yields observed with error (\(se\));

4. The variance-covariance matrix from the VAR process (\(SSZ\));

5. The intercept matrix of the physical dynamics (\(K0Z\));

6. The feedback matrix of the physical dynamics (\(K1Z\)).

The parameters \(K0Z\) and \(K1Z\) have closed-form solutions. Similarly, \(r0\) and \(se\) are derived analytically and are factored out of the log-likelihood function. In contrast, the remaining parameters, \(K1XQ\) and \(SSZ\), must be estimated numerically.

The optimization routine in MultiATSM combines the Nelder–Mead and L-BFGS-B algorithms, executed sequentially and repeated until convergence is achieved. At each iteration, the parameter vector yielding the highest likelihood is retained, enhancing robustness to local optima without resorting to full multi-start procedures. Convergence is achieved when the absolute change in the mean log-likelihood falls below a user-defined tolerance (default \(10^{-4}\)). For the bootstrap replications, the same optimization procedure is applied; however, only the Nelder–Mead algorithm is used to reduce computation time.

7 Full implementation of ATSMs

Package workflow

The complete workflow of the MultiATSM package is built around seven core functions, which together support a streamlined and modular process. An overview of these functions is provided below:

1. LabFac(): returns a list of risk factor labels used throughout the package. In particular, these labels assist in structuring sub-function inputs and generating variable and graph labels in a parsimonious manner;

2. InputsForOpt(): collects and processes the inputs needed to build the likelihood function as specified in Section 5. It estimates the model’s \(\mathbb{P}\)-dynamics and returns an object of class ATSMModelInputs, which includes print() and summary() S3 methods. The print() method summarizes model inputs and system features, while summary() reports statistics on risk factors and bond yields;

3. Optimization(): performs the estimation of the model parameters, primarily the \(\mathbb{Q}\)-dynamics, using numerical optimization. This function returns a comprehensive list of the model’s point estimates and can be computationally intensive;

4. InputsForOutputs(): an auxiliary function that compiles the necessary elements for producing numerical and graphical outputs. It also creates separate folders in the user’s Folder2Save directory to store the generated figures;

5. NumOutputs(): produces the numerical outputs as selected in Section 5.3, based on the model’s point estimates. The function returns an object of class ATSMNumOutputs, for which an autoplot() S3 method is available. This method provides a convenient way to visualize the selected graphical outputs;

6. Bootstrap(): computes confidence bounds for the numerical outputs using the bootstrap procedures defined in Section 5.3 (subsection “Bootstrap settings”). The function returns an ATSMModelBoot object, which can be accessed via the autoplot() S3 method to generate the desired graphical outputs with confidence intervals. As this step involves repeated model estimation, it may require several hours (possibly days) to complete;

7. ForecastYields(): generates bond yield forecasts and the corresponding forecast errors according to the specifications outlined in Section 5.3 (subsection “Out-of-sample forecast settings”). This function returns an object of class ATSMModelForecast, accesible via the plot() S3 method, which displays Root Mean Squared Errors (RMSEs) by country and forecast horizon.

Complete implementation

This section illustrates how to fully implement ATSMs using the MultiATSM package. A simplified two-country JPS original framework serves as the example. The implementation steps are outlined below, and a sample of graphical outputs are presented in Figures 2 – 5.

library(MultiATSM)
# 1) USER INPUTS
# A) Load database data
LoadData("CM_2024")

# B) GENERAL model inputs
ModelType <- "JPS original" 
Economies <- c("China", "Brazil") 
GlobalVar <- c("Gl_Eco_Act") 
DomVar <- c("Eco_Act") 
N <- 2  
t0_sample <- "01-05-2005" 
tF_sample <- "01-12-2019" 
OutputLabel <- "Test" 
DataFreq <-"Monthly"
Folder2Save <- NULL
StatQ <- FALSE 

# B.1) SPECIFIC model inputs
# GVAR-based models 
GVARlist <- list( VARXtype = "unconstrained", W_type = "Sample Mean", t_First_Wgvar = "2005",
                  t_Last_Wgvar = "2019", DataConnectedness = TradeFlows ) 

# JLL-based models 
JLLlist <- list(DomUnit =  "China")

# BRW inputs  
WishBC <- FALSE 
BRWlist <- within(list(Cent_Measure = "Mean", gamma = 0.05, N_iter = 250, B = 50, checkBRW = TRUE,
                       B_check = 1000, Eigen_rest = 1),  N_burn <- round(N_iter * 0.15))

# C) Decide on Settings for numerical outputs
WishFPremia <- TRUE 
FPmatLim <- c(60,120) 
                      
Horiz <- 30
DesiredGraphs <- c() 
WishGraphRiskFac <- FALSE
WishGraphYields <- FALSE
WishOrthoJLLgraphs <- FALSE

# D) Bootstrap settings
WishBootstrap <- TRUE 
BootList <- list(methodBS = 'bs', BlockLength = 4, ndraws = 5, pctg =  95)

# E) Out-of-sample forecast
WishForecast <- TRUE 
ForecastList <- list(ForHoriz = 12,  t0Sample = 1, t0Forecast = 162, ForType = "Rolling")

##########################################################################################
# NO NEED TO MAKE CHANGES FROM HERE:
# The sections below automatically process the inputs provided above, run the model 
# estimation, generate the numerical and graphical outputs, and save results.

# 2) Minor preliminary work: get the sets of factor labels
FactorLabels <- LabFac(N, DomVar, GlobalVar, Economies, ModelType)

# 3) Prepare the inputs of the likelihood function
ATSMInputs <- InputsForOpt(t0_sample, tF_sample, ModelType, Yields, GlobalMacro, 
                           DomMacro, FactorLabels, Economies, DataFreq, GVARlist, 
                           JLLlist, WishBC, BRWlist)

# 4) Optimization of the ATSM (Point Estimates)
ModelParaList <- Optimization(ATSMInputs, StatQ, DataFreq, FactorLabels, Economies, ModelType)

# 5) Numerical and graphical outputs
# a) Prepare list of inputs for graphs and numerical outputs
InputsForOutputs <- InputsForOutputs(ModelType, Horiz, DesiredGraphs, OutputLabel, StatQ, 
                                     DataFreq, WishGraphYields, WishGraphRiskFac, 
                                     WishOrthoJLLgraphs, WishFPremia, 
                                     FPmatLim, WishBootstrap, BootList, 
                                     WishForecast, ForecastList)
                                     
# b) Fit, IRF, FEVD, GIRF, GFEVD, and Term Premia
NumericalOutputs <- NumOutputs(ModelType, ModelParaList, InputsForOutputs, 
                               FactorLabels, Economies, Folder2Save)

# c) Confidence intervals (bootstrap analysis)
BootstrapAnalysis <- Bootstrap(ModelType, ModelParaList, NumericalOutputs, Economies, 
                               InputsForOutputs, FactorLabels, JLLlist, GVARlist, 
                               WishBC, BRWlist, Folder2Save)

# 6) Out-of-sample forecasting
Forecasts <- ForecastYields(ModelType, ModelParaList, InputsForOutputs, FactorLabels, 
                            Economies, JLLlist, GVARlist, WishBC, BRWlist, 
                            Folder2Save)

Figure 2: Chinese bond yield maturities with model fit comparisons. Model-fit reflects estimation using only risk-neutral (\(\mathbb{Q}\)) dynamics parameters, while Model-Implied incorporates both physical (\(\mathbb{P}\)) and risk-neutral (\(\mathbb{Q}\)) dynamics. The \(x\)-axes represent time in months and the \(y\)-axis is in natural units.

Figure 3: IRFs from the Brazilian bond yields to global economic activity. Size of the shock is one-standard deviation. The black lines are the point estimates. Gray dashed lines are the bounds of the 95% confidence intervals and the green lines correspond to the median of these intervals. The \(x\)-axes are expressed in months and the \(y\)-axis is in natural units.

Figure 4: FEVD from the Brazilian bond yield with maturity 60 months. The \(x\)-axis represents the forecast horizon in months and the \(y\)-axis is in natural units.

Figure 5: Chinese sovereign yield curve decomposition showing (i) expected future short rates and (ii) term premia components. The \(x\)-axis represents time in months and the \(y\)-axis is expressed in percentage points.

8 Concluding remarks

The MultiATSM package aims to advance yield curve (term structure) modelling within the R programming environment. It provides a comprehensive yet user-friendly toolkit for practitioners, academics, and policymakers, featuring estimation routines and generating detailed outputs across several macrofinance model classes. This allows for an in-depth exploration of the relationship between the real economy developments and the fixed income markets.

The package covers eight classes of macrofinance term structure models, all built upon the single-country unspanned macroeconomic risk framework of Joslin et al. (2014), which is also extended to a multicountry setting. Additional multicountry variants based on Jotikasthira et al. (2015) and Candelon and Moura (2024) are included, incorporating, respectively, a dominant economy and a GVAR structure to model cross-country interdependence.

Each model class provides analytical outputs that offer insight into term structure dynamics, including plots of model fit, risk premia, impulse responses, and forecast error variance decompositions.The MultiATSM package also offers bootstrap procedures for confidence interval construction and out-of-sample forecasting of bond yields.

Acknowledgments

I thank the editor, Rob Hyndman, and an anonymous referee for several helpful comments. I am also grateful to Bertrand Candelon, Adhir Dhoble and Gustavo Torregrosa for many insightful discussions. An earlier version of this paper circulated under the title MultiATSM: An R Package for Arbitrage-Free Multicountry Affine Term Structure of Interest Rate Models with Unspanned Macroeconomic Risk and was part of the author’s PhD dissertation at UCLouvain (Moura 2022). The views expressed in this paper are those of the author and do not necessarily reflect those of Banco de Mexico.

9 Appendix

A: Supplementary functions

Importing data from Excel files

The MultiATSM package also provides an automated procedure for importing data from Excel files via Load_Excel_Data() and preparing the risk factor database used directly in the model estimation. To ensure compatibility with the package functions, the following requirements must be met:

1. Databases must be organized in separate Excel files: one for unspanned factors and another for term structure data. For GVAR-based models, a third file containing the interdependence measures is also required;
2. Each Excel file should include one tab per country. In the case of unspanned factors, an additional tab must be included for the global variables if the user opts to incorporate them;
3. Variable names must be identical across all tabs within each file.

An example Excel file meeting these requirements is provided with the package. Below is an example of how to import the data from excel and construct the input list to be supplied:

MacroData  <- Load_Excel_Data(system.file("extdata", "MacroData.xlsx", 
                                          package = "MultiATSM"))
YieldsData <- Load_Excel_Data(system.file("extdata", "YieldsData.xlsx", 
                                          package = "MultiATSM"))

ModelType <- "JPS original"
Initial_Date <- "2006-09-01"
Final_Date <- "2019-01-01"
DataFrequency <- "Monthly"
GlobalVar <- c("GBC", "VIX")
DomVar <- c("Eco_Act", "Inflation", "Com_Prices", "Exc_Rates")
N <- 3
Economies <- c("China", "Mexico", "Uruguay", "Brazil", "Russia")

These inputs are used to construct the RiskFactorsSet variable, which holds the full collection of risk factors required by the model.

FactorLabels <- LabFac(N, DomVar, GlobalVar, Economies, ModelType)
RiskFactorsSet <- DataForEstimation(Initial_Date, Final_Date, Economies, N, FactorLabels,
                                    ModelType, DataFrequency, MacroData, YieldsData)

Transition matrix and star factors

To construct the transition matrix for GVAR specifications, the user can employ Transition_Matrix(). This function requires:

1. Data selection: choose proxies for cross-country interdependence.

2. Time frame: specify the sample’s start and end dates.

3. Dependence measure: select from:

• Time-varying (dynamic weights)
• Sample Mean (static average)
• A numeric scalar (fixed-year snapshot).

data("TradeFlows")
t_First <- "2006"
t_Last <- "2019"
Economies <- c("China", "Brazil", "Mexico", "Uruguay")
type <- "Sample Mean"
W_gvar <- Transition_Matrix(t_First, t_Last, Economies, type, TradeFlows)

Note that if data is missing for any country in a given year, the corresponding transition matrix will contain only NAs.

A more flexible approach to modelling interdependence is to allow the transition matrix to vary over time. In this case, the star factors are constructed using trade flow weights specific to each year, adjusting the corresponding year’s risk factors accordingly. To enable this feature, users must set the type argument to Time-varying and specify the same year for both the initial and final periods in the transition matrix. This indicates that the trade weights from that particular year is used when solving the GVAR system (i.e., in the construction of the link matrices, see Equation (11)).

B: Additional theoretical considerations

Bond yield decomposition

The MultiATSM package allows for the calculation of two risk compensation measures: term premia and forward premia. Assume that an \(n\)-maturity bond yield can be decomposed into two components: the expected short-rate (\(\mathrm{Exp}_{i,t}^{(n)}\)) and term premia (\(\mathrm{TP}_{i,t}^{(n)}\)). Technically: \[ y_{i,t}^{(n)} = \mathrm{Exp}_{i,t}^{(n)} + \mathrm{TP}_{i,t}^{(n)} \text{.} \] In the package’s standard form, the expected short rate term is computed from time \(t\) to \(t+n\), which represents the bond’s maturity: \(\mathrm{Exp}_{i,t}^{(n)} = \sum_{h=0}^{n} E_t[y_{i, t+h}^{(1)}]\). Alternatively, the decomposition for the forward rates (\(f_{i,t}^{(n)}\)) is \(f_{i,t}^{(n)} = \sum_{h=m}^{n} E_t[y_{i,t+h}^{(1)}] + \mathrm{FP}_{i,t}^{(n)}\) where \(\mathrm{FP}_{i,t}^{(n)}\) corresponds to the forward premia. In this case, the user must specify TRUE if the computation of forward premia is desired, or FALSE otherwise. If set to TRUE, the user must also provide a two-element numerical vector containing the maturities corresponding to the starting and ending dates of the bond maturity. Example:

    WishFPremia <- TRUE  
    FPmatLim <- c(60, 120)

C: Replication of existing research

Joslin, Priebisch and Singleton (2014)

The dataset used in this replication was constructed by Bauer and Rudebusch (2017) (henceforth BR, 2017) and is available on Bauer’s website. In their paper, BR (2017) investigate whether macrofinance term structure models are better suited to the unspanned macro risk framework of JPS (2014) or to earlier, traditional spanned settings such as Ang and Piazzesi (2003). To that end, BR (2017) replicate selected empirical results from JPS (2014). The corresponding R code is also available on Bauer’s website.

Using the dataset from BR (2017), the code below applies the MultiATSM package to estimate the key ATSM parameters following the JPS original modelling setup.

# 1) INPUTS
# A) Load database data
LoadData("BR_2017")

# B) GENERAL model inputs
ModelType <- "JPS original"

Economies <- c("US") 
GlobalVar <- c() 
DomVar <- c("GRO", "INF")
N <- 3 
t0_sample <- "January-1985"
tF_sample <- "December-2007"
DataFreq <- "Monthly" 
StatQ <- FALSE 

# 2) Minor preliminary work
FactorLabels <- LabFac(N, DomVar, GlobalVar, Economies, ModelType) 
Yields <- t(BR_jps_out$Y)
DomesticMacroVar <- t(BR_jps_out$M.o)
GlobalMacroVar <- c()

# 3) Prepare the inputs of the likelihood function
ATSMInputs <- InputsForOpt(t0_sample, tF_sample, ModelType, Yields, GlobalMacroVar, 
                           DomesticMacroVar, FactorLabels, Economies, DataFreq)

# 4) Optimization of the model
ModelPara <- Optimization(ATSMInputs, StatQ, DataFreq, FactorLabels, Economies, ModelType)

The tables below compare the ATSM parameter estimates generated from BR (2017) and the MultiATSM. Table 2 reports the risk-neutral parameters. While the values presented do not match exactly, the differences are well within convergence tolerance and are arguably economically negligible. Table 3, by contrast, contains parameters related to the model’s time-series dynamics. As these are derived in closed form, the estimates are exactly the same under both specifications.

Table 2: \(Q\)-dynamics parameters

	MultiATSM	BR (2017)
\(r_0\)	0.0006	−0.0002
\(\lambda_1\)	0.9967	0.9968
\(\lambda_2\)	0.9149	0.9594
\(\lambda_3\)	0.9149	0.8717
Note:
λ’s are the eigenvalues from the risk-neutral feedback matrix and r₀ is the long-run mean of the short rate under Q.

Table 3: \(P\)-dynamics parameters

	K0Z	K1Z
		PC1	PC2	PC3	GRO	INF
BR (2017)
PC1	0.0781	0.9369	−0.0131	−0.0218	0.1046	0.1003
PC2	0.0210	0.0058	0.9781	0.1703	−0.1672	−0.0402
PC3	0.1005	−0.0104	−0.0062	0.7835	−0.0399	0.0437
GRO	0.0690	−0.0048	0.0180	−0.1112	0.8818	−0.0025
INF	0.0500	0.0018	0.0064	−0.0592	0.0277	0.9859
MultiATSM
PC1	0.0781	0.9369	−0.0131	−0.0218	0.1046	0.1003
PC2	0.0210	0.0058	0.9781	0.1703	−0.1672	−0.0402
PC3	0.1005	−0.0104	−0.0062	0.7835	−0.0399	0.0437
GRO	0.0690	−0.0048	0.0180	−0.1112	0.8818	−0.0025
INF	0.0500	0.0018	0.0064	−0.0592	0.0277	0.9859
Note:
\(K0Z\) is the intercept and \(K1Z\) is the feedback matrix from the \(P\)-dynamics.

For replicability, it is important to note that the physical dynamics results reported in Table 3 using MultiATSM rely on the principal component weights provided by BR (2017). Such a matrix is simply a scaled-up version of the one provided by the function pca_weights_one_country() of the current package. Accordingly, despite the numerical differences on the weight matrices, both methods generate time series of spanned factors which are perfectly correlated. Another difference between the two approaches relates to the construction of the log-likelihood function: while in the BR (2017) code this is expressed in terms of a portfolio of yields, the MultiATSM package generates this same input directly as a function of observed yields (i.e. both procedures lead to equivalent log-likelihood vales up to the Jacobian term).

Additionally, it is worth highlighting that the standard deviations for the portfolios of yields observed with errors are nearly identical, matching to seven decimal places: 0.0000546 for MultiATSM and 0.0000550 for BR (2017).

Candelon and Moura (2024)

The multicountry framework introduced in Candelon and Moura (2024) enhances the tractability of large-scale ATSMs and deepens our understanding of the global economic mechanisms driving domestic yield curve fluctuations. This framework also generates more precise model estimates and enhances the forecasting capabilities of these models. This novel setup, embodied by the GVAR multi model class, is benchmarked against the findings of Jotikasthira et al. (2015), which are captured by the JLL original model class. The paper showcases an empirical illustration involving China, Brazil, Mexico, and Uruguay.

# 1) INPUTS
# A) Load database data
LoadData("CM_2024")

# B) GENERAL model inputs
ModelType <- "GVAR multi" 
Economies <- c("China", "Brazil", "Mexico", "Uruguay")
GlobalVar <- c("Gl_Eco_Act", "Gl_Inflation")
DomVar <- c("Eco_Act", "Inflation")
N <- 3
t0_sample <- "01-06-2004"
tF_sample <- "01-01-2020"
OutputLabel <- "CM_jfec"
DataFreq <-"Monthly"
StatQ <- FALSE

# B.1) SPECIFIC model inputs
# GVAR-based models 
GVARlist <- list( VARXtype = "unconstrained", W_type = "Sample Mean", t_First_Wgvar = "2004",
                  t_Last_Wgvar = "2019", DataConnectedness = TradeFlows )

# JLL-based models 
JLLlist <- list(DomUnit =  "China")

# BRW inputs
WishBC <- TRUE
BRWlist <- within(list(Cent_Measure = "Mean", gamma = 0.001, N_iter = 200, B = 50, checkBRW = TRUE,
                       B_check = 1000, Eigen_rest = 1),  N_burn <- round(N_iter * 0.15))

# C) Decide on Settings for numerical outputs
WishFPremia <- TRUE
FPmatLim <- c(24,36)

Horiz <- 25
DesiredGraphs <- c("GIRF", "GFEVD", "TermPremia")
WishGraphRiskFac <- FALSE
WishGraphYields <- TRUE
WishOrthoJLLgraphs <- TRUE

# D) Bootstrap settings
WishBootstrap <- FALSE 
BootList <- list(methodBS = 'bs', BlockLength = 4, ndraws = 1000, pctg =  95)

# E) Out-of-sample forecast
WishForecast <- TRUE
ForecastList <- list(ForHoriz = 12,  t0Sample = 1, t0Forecast = 100, ForType = "Rolling")

# 2) Minor preliminary work: get the sets of factor labels and  a vector of common maturities
FactorLabels <- LabFac(N, DomVar, GlobalVar, Economies, ModelType)

# 3) Prepare the inputs of the likelihood function
ATSMInputs <- InputsForOpt(t0_sample, tF_sample, ModelType, Yields, GlobalMacro, 
                           DomMacro, FactorLabels, Economies, DataFreq, 
                           GVARlist, JLLlist, WishBC, BRWlist)

# 4) Optimization of the ATSM (Point Estimates)
ModelParaList <- Optimization(ATSMInputs, StatQ, DataFreq, FactorLabels, Economies, ModelType)

# 5) Numerical and graphical outputs
# a) Prepare list of inputs for graphs and numerical outputs
InputsForOutputs <- InputsForOutputs(ModelType, Horiz, DesiredGraphs, OutputLabel, StatQ, 
                                     DataFreq, WishGraphYields, WishGraphRiskFac, 
                                     WishOrthoJLLgraphs, WishFPremia, FPmatLim, 
                                     WishBootstrap, BootList, WishForecast, 
                                     ForecastList)
                                    
# b) Fit, IRF, FEVD, GIRF, GFEVD, and Term Premia
NumericalOutputs <- NumOutputs(ModelType, ModelParaList, InputsForOutputs, 
                               FactorLabels, Economies)

# c) Confidence intervals (bootstrap analysis)
BootstrapAnalysis <- Bootstrap(ModelType, ModelParaList, NumericalOutputs, Economies,
                               InputsForOutputs, FactorLabels, JLLlist, GVARlist, 
                               WishBC, BRWlist)
                               
# 6) Out-of-sample forecasting
Forecasts <- ForecastYields(ModelType, ModelParaList, InputsForOutputs, FactorLabels, 
                            Economies, JLLlist, GVARlist, WishBC, BRWlist)

Candelon and Moura (2023)

In this paper, Candelon and Moura (2023) investigate the underlying factors that shape the sovereign yield curves of Brazil, India, Mexico, and Russia during the COVID\(-19\) pandemic crisis. The study adopts a GVAR multi approach to capture the complex global macrofinancial, and especially health-related interdependencies during the latest pandemic.

# 1) INPUTS
# A) Load database data
LoadData("CM_2023")

# B) GENERAL model inputs
ModelType <- "GVAR multi"
Economies <- c("Brazil", "India", "Russia", "Mexico")
GlobalVar <- c("US_Output_growth", "China_Output_growth", "SP500")
DomVar <- c("Inflation","Output_growth", "CDS", "COVID")
N <- 2
t0_sample <- "22-03-2020"
tF_sample <- "26-09-2021"
OutputLabel <- "CM_EM"
DataFreq <-"Weekly"
StatQ <- FALSE

# B.1) SPECIFIC model inputs
# GVAR-based models 
GVARlist <- list(VARXtype = "constrained: COVID", W_type = "Sample Mean", 
                 t_First_Wgvar = "2015", t_Last_Wgvar = "2020", 
                 DataConnectedness = TradeFlows_covid)

# BRW inputs  
WishBC <- FALSE

# C) Decide on Settings for numerical outputs
WishFPremia <- TRUE
FPmatLim <- c(47,48)

Horiz <- 12
DesiredGraphs <- c("GIRF", "GFEVD", "TermPremia")
WishGraphRiskFac <- FALSE
WishGraphYields <- TRUE
WishOrthoJLLgraphs <- FALSE

# D) Bootstrap settings
WishBootstrap <- TRUE 
BootList <- list(methodBS = 'bs', BlockLength = 4, ndraws = 100, pctg =  95)

# 2) Minor preliminary work: get the sets of factor labels and  a vector of common maturities
FactorLabels <- LabFac(N, DomVar, GlobalVar, Economies, ModelType)

# 3) Prepare the inputs of the likelihood function
ATSMInputs <- InputsForOpt(t0_sample, tF_sample, ModelType, Yields_covid, GlobalMacro_covid, 
                           DomMacro_covid, FactorLabels, Economies, DataFreq, GVARlist) 
                           
# 4) Optimization of the ATSM (Point Estimates)
ModelParaList <- Optimization(ATSMInputs, StatQ, DataFreq, FactorLabels, Economies, ModelType)

# 5) Numerical and graphical outputs
# a) Prepare list of inputs for graphs and numerical outputs
InputsForOutputs <- InputsForOutputs(ModelType, Horiz, DesiredGraphs, OutputLabel, StatQ, 
                                     DataFreq, WishGraphYields, WishGraphRiskFac, 
                                     WishOrthoJLLgraphs, WishFPremia, FPmatLim,
                                     WishBootstrap, BootList)

# b) Fit, IRF, FEVD, GIRF, GFEVD, and Term Premia
NumericalOutputs <- NumOutputs(ModelType, ModelParaList, InputsForOutputs, FactorLabels, 
                               Economies)

# c) Confidence intervals (bootstrap analysis)
BootstrapAnalysis <- Bootstrap(ModelType, ModelParaList, NumericalOutputs, Economies, 
                               InputsForOutputs, FactorLabels, 
                               JLLlist = NULL, GVARlist)

10 Supplementary materials

Supplementary materials are available in addition to this article. It can be downloaded at RJ-2025-044.zip

11 CRAN packages used

MultiATSM, YieldCurve, fBonds, statespacer, simStateSpace, vars, MTS, svars, bigtime, BigVAR, Spillover, BGVAR

12 CRAN Task Views implied by cited packages

Econometrics, Finance, TimeSeries

---

1. Specifically, the referred loadings are \(a_{i,n+1}(\Theta _{n+1}) = a_{i,n}(\Theta _{n}) + b_{i,n}(\Theta _{n}) \mu^{Q}_{i,X} + \frac{1}{2} b_{i,n}(\Theta _{n}) \Gamma_{i,X} \Gamma_{i,X}' b_{i,n}(\Theta _{n})' - \delta_{i,0}\) and \(b_{i, n+1}=b_{i,n}\Phi^{Q}_{i,X} - \delta_{i,1}\), considering that the boundary conditions are \(a_{i,1}(\Theta _1)=-\delta_{i,0}\) and \(b_{i,1}(\Theta_1)=-\delta_{i,1}\). These expressions assume that the Radon–Nikodym derivative, which maps the risk-neutral measure to the physical measure, follows a log-normal process, and that the market price of risk is time-varying and affine in \(X_t\). See Ang and Piazzesi (2003) for a detailed derivation of these expressions.↩︎

2. Noticeably, in the context of the MultiATSM package, the model type JLL No DomUnit is the only exception (see Section 3).↩︎

3. It is worth emphasizing that, although the DataForEstimation() and InputsForOpt() functions in the package accept inputs with differing maturities, their outputs are standardized to a common set of yields.↩︎