---
title: "MultiATSM package - General Guidelines"
author: Rubens Moura
date: "`r Sys.Date()`"
output: 
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
    number_sections: true
    df_print: paged
    toc: yes
    toc_depth: 2
    fig_caption: yes
vignette: >
  %\VignetteIndexEntry{Vignette for package MultiATSM}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc} 
bibliography: references.bib
---

This document aims at providing general guidance on the use of the *MultiATSM* package available at the [CRAN repository](https://cran.r-project.org/). 

```{r}
library(MultiATSM)
```

The *MultiATSM* package provides several model outputs from various single and multi-country affine term structure of interest models (ATSMs). All the frameworks of this package are based on the unspanned economic risk framework from @JoslinPriebschSingleton2014 (JPS, 2014). In essence, these models assume the absence of arbitrage opportunities, consider a linear state space representation of the yield curve dynamics and offer a tractable approach to simultaneously combine the traditional yield curve factors (spanned factors) along with economic and financial variables (unspanned factors). 

Due to the peculiar features of the JPS specification, an efficient estimation of the parameters governing the risk-neutral ($Q$-measure) and the physical ($P$-measure) probability measures can be carried out rather independently. The only exception is the variance-covariance matrix (sigma) term which is a common element to both the $P$ and the $Q$ density functions.

In its current version, the *MultiATSM* package can generate the outputs of 8 different classes of ATSMs. More specifically, the *MultiATSM* includes models in which the estimation is done on a *country-by-country* basis - as in JPS (2014)  - or *jointly* for all the countries of the economic system - as in @JotikasthiraLeLundblad2015 (JLL, 2015) and @CandelonMoura2021 (CM, 2021).  In all cases, the risk factor dynamics under the $P$-measure follow a VAR(1) model of some sort. The table below summarizes the general features of each one of the models available at this package.   

```{r ModFea, message=FALSE, echo=FALSE}
ModelLabels <- c("JPS", "JPS jointP", "VAR jointQ", "GVAR sepQ", "GVAR jointQ", 
                 "JLL original", "JLL NoDomUnit", "JLL jointSigma")

# Rows
Tab <- data.frame(matrix(nrow = length(ModelLabels), ncol = 0)) 
rownames(Tab) <- ModelLabels

# Empty columns
EmptyCol <- c("", "", "", "", "", "", "", "") 
Tab$EmptyCol0 <- EmptyCol
# P-dynamics + 2 empty spaces
Tab$PdynIndUnco <- c("x", "", "", "", "", "", "", "")
Tab$PdynIndCo <- c("", "", "", "", "", "", "", "")
Tab$PdynJointUnco <- c("", "x", "x", "", "", "", "", "")
Tab$PdynJointJLL <- c("", "", "", "", "", "x", "x", "x")
Tab$PdynJointGVAR <- c("", "", "", "x", "x", "", "", "")
Tab$EmptyCol1 <- EmptyCol
Tab$EmptyCol2 <- EmptyCol
# Q-dynamics + 2 empty spaces
Tab$QdynInd <- c("x", "x", "", "x", "", "", "", "")   
Tab$QdynJoint <- c("", "", "x", "", "x", "x", "x", "x") 
Tab$EmptyCol3 <- EmptyCol
Tab$EmptyCol4 <- EmptyCol
# Sigma + 2 empty spaces
Tab$Ponly <-  c("", "", "", "", "", "x", "x", "")
Tab$PandQ <- c("x", "x", "x", "x", "x", "", "", "x")
Tab$EmptyCol5 <- EmptyCol
Tab$EmptyCol6 <- EmptyCol
# Dominant Unit
Tab$DomUnit <- c("", "", "", "", "", "x", "", "x")

# Adjust column names
ColNames <- c("","","","","JLL", "GVAR", "", "", "", "", "", "", "","", "", "","")
colnames(Tab) <- ColNames


# Generate the table
suppressWarnings(library(magrittr))


kableExtra::kbl(Tab, align = "c", caption = "Model Features") %>%
  kableExtra::kable_classic("striped", full_width = F)  %>%
  kableExtra::row_spec(0, font_size = 10) %>%
  kableExtra::add_header_above(c(" "=2, "Unrestricted" = 1, "Restricted" = 1, "Unrestricted" = 1, "Restricted" = 2, " " = 11)) %>%
 kableExtra::add_header_above(c(" "=2, "Individual" = 2, "Joint" = 3, " "=2, "Individual" = 1, "Joint" = 1, " "=2, "P only" = 1, "P and Q" = 1, " " = 3))  %>%
  kableExtra::add_header_above(c( " "=2, "P-dynamics"= 5, " "=2, "Q-dynamics"= 2, " "=2, "Sigma matrix estimation" = 2, " "=2, "Dominant Country"=1), bold = T) %>%
kableExtra::pack_rows("Unrestricted VAR", 1, 3 , label_row_css = "background-color: #666; color: #fff;")  %>%
kableExtra::pack_rows("Restricted VAR (GVAR)", 4, 5, label_row_css = "background-color: #666; color: #fff;") %>%
kableExtra::pack_rows("Restricted VAR (JLL)", 6, 8, label_row_css = "background-color: #666; color: #fff;")
```

Some aspects of the multi-country frameworks are worth highlighting. As for the models based on the setup of JLL (2015), the version "JLL original" follows closely the seminal work of JLL (2015), i.e., it is assumed an economic cohort containing a worldwide dominant economy and a set of smaller countries, in addition to the estimation of the sigma matrix be performed exclusively under the $P$-measure. The two other alternative versions assume the absence of a dominant country ("JLL NoDomUnit") and the estimation of sigma under both the $P$ and $Q$ measures ("JLL jointSigma"), as in the standard JPS (2014) model. As for the remaining multi-country ATSMs, it is considered that the dynamics of the risk factors under the $P$-measure evolves according to a GVAR model. The version labeled "GVAR jointQ" is the one presented in CM (2021).

In what follows, the focus will be devoted to describe the various pieces that are necessary to implement ATSMs. Specifically, Section \@ref(S:SectionData) describes the data-set available at this package and the set of functions that are useful to retrieve data from Excel files. Section \@ref(S:SectionInputs) exposes the necessary inputs that have to be specified by the user. In Section \@ref(S:SectionEstimation) the estimation procedure is detailed. Section \@ref(S:SectionReplication) shows how to use the *MultiATSM* package to estimate ATSMs from scratch as well as to closely reproduce some empirical features of academic papers. To ease the exposition, throughout the next sections, the database used by CM (2021) is employed to illustrate the use of the various functions available in this package.

# Data {#S:SectionData}
 
## Package data-set{#S:SectionDataPackage}
The *MultiATSM* package contains the four data-sets used in CM (2021). The first set of data comprises several time series of zero-coupon bond yields from four emerging markets: China, Brazil, Mexico, and Uruguay. It is worth noting that this package requires *(i)* for estimation purposes, the maturities of the bond yields to be the same for all countries (although the function *DataForEstimation* handles different maturities across countries as inputs, the outputs generated by this same function are bond yields of the same maturities for all the economies); *(ii)* bond yields to be expressed in percentage terms (and not in basis points) per annum. The *MultiATSM* package does not support a routine to bootstrap zero-coupon yields from coupon bonds. This data manipulation procedure must be handled by the user herself, if necessary.  
```{r}
data('CM_Yields')
```

The second group of data concerns the time series of the risk factors. Specifically, along the terminology defined in JPS (2014), this data-set contains *(i)* country-specific spanned factors (denoted by level, slope, and curvature - see Section \@ref(S:SpaFac) for a more detailed description) and *(ii)* an array of country-specific and global unspanned factors (namely, some measure of economic growth and inflation). Similarly to the case of the bond yield data, the measures of economic and financial variables must be constructed by the user.
```{r}
data('CM_Factors')
```

The two last blocks of data are necessary for the estimation of the GVAR-based models. The trade flows database
presents the sum of the value of all goods imports and exports between any two countries of the sample on a yearly basis since 1948. All the values are free on board and are expressed in U.S. dollars. These data are used to construct the transition matrix of the GVARs models.  
```{r}
data('CM_Trade')
```

The GVAR factors database casts country-specific lists of all the factors that are used in the estimation of each country's VARX(1,1,1). A specific function for computing the star risk factors is detailed in the Section \@ref(S:SectionGVAR). 
```{r}
data('CM_Factors_GVAR')
```

## Import data from Excel files {#S:SectionImport}
This package also offers an automatized procedure to extract data from Excel files and to, subsequentially, prepare the risk factor databases that are directly used in the model estimation. The use of the package functions requires that the databases *(i)* are constructed in separate Excel files for the unspanned factors, the term structure data, and the measures of interdependence (for the GVAR-based models); *(ii)* contain, in each Excel file, one separate tab per country (in addition to the global variables, in the case of the unspanned factors' database); and *(iii)* have identical variable labels across all the tabs within each file. For illustration, see the Excel file available at the package.
One example of list of inputs to be provided is, for instance,
```{r}
Initial_Date <- "2006-09-01" # Format "yyyy-mm-dd"
Final_Date <- "2019-01-01" # Format "yyyy-mm-dd"
DataFrequency <- "Monthly"
GlobalVar <- c("GBC", "VIX") # Global Variables
DomVar <- c("Eco_Act", "Inflation", "Com_Prices", "Exc_Rates") #  Domestic variables
N <-  3 # Number of spanned factors per country
Economies <- c("China", "Mexico", "Uruguay", "Brazil", "Russia")
ModelType <- "JPS"
```

Based on these inputs, one can construct the variable *ZZfull* which contains the complete set of the model risk factors.
```{r}
FactorLabels <- LabFac(N, DomVar, GlobalVar, Economies, ModelType)
ZZfull <-DataForEstimation(Initial_Date, Final_Date, Economies, N, FactorLabels, ModelType, DataFrequency)
```

# Required user inputs {#S:SectionInputs}

## Basic inputs used in the model estimation 

In order to estimate any ATSM of this package, the user needs to specify several general model inputs, namely:

- **Model type**: a string-vector containing the label of the model to be estimated (as per described in Table \@ref(tab:ModFea)); 
- **Frequency of the data**: a string-vector that specifies the frequency of the time series data. The available options are: *"Annually"*, *"Quarterly"*, *"Monthly"*, *"Weekly"*, *"Daily Business Days"*, and *"Daily All Days"*;
- **Economies**: a string-vector containing the names of the economies which are part of the economic system;
- **Number of spanned factors (N)**: number of country-specific spanned factors. In this package, *N* is assumed to be equal across countries;
- **Global variables**: a string-vector containing the labels of the global unspanned factors used in the model estimation;
- **Domestic variables**: a string-vector containing the labels of the domestic unspanned factors used in the model estimation;
- **Bond yield set of maturities**: a numerical-vector containing the maturities of the yields used in the model estimation. The units of this vector must be expressed in *years*. Alternatively, the user can rely on the *Maturities* routine to extract the bond maturities from the time series of bond yields (see the documentation of the *Maturities* function). In this case, the user must specify the **unit in which the label of the bond yields** are expressed, e.g. if the label of the 10-years bond is selected as *Y120M* (*Y10y*), then the option *Month* (*Year*) must be selected;
- **Factor labels**: a string-list based which contains the labels of all the variables present in the model. This input is particularly useful for constructing the plots. The routine *LabFac* provides the desired list of factor labels;
- **Stationarity constraint under the $Q$-dynamics**: a binary variable which takes value equal to "1" if the user wishes the largest eigenvalue under the $Q$-measure to be strictly smaller than 1. Otherwise, this same variable must be set to "0";
- **Output label**: a single element string-vector which contains the of the name that appears in the file name that stores the model outputs.


One possible example of the basic model inputs is
```{r}
ModelType <- "JPS"
DataFrequency <- "Monthly"
Economies <- c("China", "Brazil", "Mexico", "Uruguay")
N <- 3 # number of spanned factors per country
GlobalVar <- c("GBC", "CPI_OECD") # Global variables
DomVar <- c("Eco_Act", "Inflation") # Domestic variables 
mat <- c(0.25, 0.5, 1, 3, 5, 10) # vector of maturities
FactorLabels <- LabFac(N, DomVar,GlobalVar, Economies, ModelType) 
StationarityUnderQ <- 0 # 1 = set stationarity condition under  the Q; 0 = no stationarity condition
OutputLabel <- "Model_demo" # output label
```


### GVARinputs and JLLinputs {#S:SectionGVARJLLinputs}
Some additional inputs are required if the user intends to estimate the GVAR-ATSM or JLL-ATSM related setups. For a convenient use, the extra outputs should be stored in a separate list for each model. The GVAR model list of inputs must contain:

- **Economies**: a string-vector containing the names of the economies as described above;  
- **GVAR list of risk factors**:  a list of risk factors sorted by country in addition to the global variables. See function *DatabasePrep*;
- **VARX type**: a string-vector containing the desired estimated form of the VARX(1,1,1). Two possibilities are available. The *"unconstrained"* form estimates the model without any constrains, using standard OLS regressions for each model equation. The *"constrained: Spanned Factors"* form prevents foreign-spanned-factors to impact any domestic factor in the feedback matrix, whereas *"constrained: '* followed by the name of the risk factor restricts this same factor to be influenced only by its own lagged values and the lagged values of its own star variables. In the last two cases, the VARX(1,1,1) is estimated by restricted least squares.


```{r}
GVARinputs <- list()
GVARinputs$Economies <- Economies
GVARinputs$GVARFactors <- FactorsGVAR
GVARinputs$VARXtype <- "constrained: Inflation"
```

Furthermore, the estimation of the GVAR requires the specification of the necessary inputs to build the **transition matrix** (see Section \@ref(S:SectionGVAR) for more details). One example of this list of inputs is
```{r}
t_First <- "2000" # First year of the sample
t_Last <-  "2015" # Last year of the sample
W_type <- 'Sample Mean' 
```


Concerning the JLL-ATSM frameworks, the list of inputs includes:  

- **JLL type**: a string-vector containing the label of the JLL model to be estimated (as described in Table \@ref(tab:ModFea));
- **Economies**: a string-vector containing the names of the economies as described above;
- **Dominant Country**: a string-vector containing the name of the economy which is assigned as the dominant country (applicable for the *JLL original* and *JLL jointSigma* models) or "None" (applicable for the *JLL NoDomUnit*).
 - **Wish the estimation of the sigmas matrices**: this is a binary variable which assumes value equal to 1 if the user wishes the estimation of all JLL sigma matrices (i.e. variance-covariance and the Cholesky factorization matrices) and, 0 otherwise. The estimation of the sigma matrices can take several minutes;
- **Sigma of the non-orthogonalized variance-covariance matrix**: to save time, the user may provide the variance-covariance matrix from the non-orthogonalized dynamics. Otherwise, this input should be set as *NULL*.

One possible list of JLL inputs is as follows:
```{r}
  JLLinputs <- list()
  ModelType <- "JLL original"  
  JLLinputs$Economies <- Economies
  JLLinputs$DomUnit <- "China"
  JLLinputs$WishSigmas <- 1 
  JLLinputs$SigmaNonOrtho <- NULL
  JLLinputs$JLLModelType <- ModelType
```


## Additional inputs for numerical and graphical outputs
Once the model parameters from the ATSM have been estimated, the *MultiATSM* package enables the numerical and graphical compilation of the following additional outputs:

- Model fit of the bond yields;
- Impulse response functions (IRFs);
- Forecast error variance decompositions (FEVDs);
- Generalized impulse response functions (GIRFs);
- Generalized forecast error variance decompositions (GFEVDs).

For the IRFs, GIRFs, FEVDs or GFEVDs, a horizon of analysis has to be specified, e.g.:
```{r}
Horiz <- 100
```

For the graphical outputs, the user must indicate the desired types of graphs in a string-based vector:
```{r}
DesiredGraphs <- c("Fit", "GIRF", "GFEVD") # Available options are: "Fit", "IRF", "FEVD", "GIRF", "GFEVD"
```

Moreover, the user must select the types of variables of interest (yields, risk factors or both) and, for the JLL type of models, whether the orthogonalized version should be additionally included   
```{r}
WishGraphRiskFac <- 0 #   YES: 1; No = 0.
WishGraphYields <- 1 #    YES: 1; No = 0.
WishOrthoJLLgraphs <- 0 # YES: 1; No = 0.
```

The function *InputsForOutputs* can provide some guidance for customizing the features of the wished outputs. Conditional to these settings, individual folders are created at your *temporary directory* to store the different types of the desired graphical outputs. 

### Bootstrap settings
If the user intends to generate confidence bands through some bootstrap procedure, an additional list of inputs is required. Specifically:

- **methodBS**: the desired bootstrap procedure. Available options are *(i)* standard residual bootstrap ('bs'); *(ii)* wild bootstrap ('wild') and block bootstrap ('block'). If the latest is selected, then the block length must be indicated;
- **ndraws**:  the number of bootstrap draws;
- **pctg**: the confidence level expressed in percentage points.

```{r}
Bootlist <- list()
Bootlist$methodBS <- 'block' 
Bootlist$BlockLength <- 4 
Bootlist$ndraws <-  3  
Bootlist$pctg   <-  95 
```


### Out-of-sample forecast settings
To perform out-of-sample forecasts of bond yields, the following list-based features have to be detailed:

- **ForHoriz**: number of periods-ahead that the forecasts are to be generated;
- **t0Sample**: time-dimension index of the first observations belonging to the information set;
- **t0Forecast**: time-dimension index of the last observation of the information set used to perform the first forecast set.

```{r}
ForecastList <- list()
ForecastList$ForHoriz <- 12 # forecast horizon
ForecastList$t0Sample <- 1 # initial sample date
ForecastList$t0Forecast <- 70 # last sample date for the first forecast
```



# Model Estimation {#S:SectionEstimation}
Having gathered bond yields and other economic time series data, the ATSM estimation is carried out in three steps. Firstly, one needs to provide the number of country-specific spanned factors to be included in the global ATSM. Secondly, it is decided on the form of the risk factor dynamics under the $P$-measure. Finally, the user needs to select general features for the model optimization.

## Spanned Factors {#S:SpaFac}
The spanned factors are yield-related variables that are used to fit the cross-section dimensions of the term structures. Typically, the spanned factors of country $i$, $P_{i,t}$, are computed as the *N* first principal components (PCs) of the set of observed bond yields. Formally, $P_{i,t}$ is constructed as $P_{i,t} = w_i Y_{i,t}$ where $w_i$ is the PC weight matrix and $Y_{i,t}$ is a country-specific column-vector of yields with increasing maturities. 

For *N=3*, the spanned factors are traditionally interpreted as level, slope, and curvature. The nature of such interpretability results from the features of the PC weight matrix as illustrated below: 

```{r}
w <- pca_weights_one_country(Yields, Economy = "Uruguay") 
```

In matrix *w*, each row stores the weights used for constructing each spanned factor. The entries of the first row are linked to the composition of the *level* factor in that they load rather equally on all yields. Accordingly, high (low) values of the level factor indicate an overall high (low) value of yields across all maturities. In the second row, the weights monotonically increase with the maturities and, therefore, they capture the degree of steepness (*slope*) of the term structure. High slope factor values imply a steep yield curve, whereas low ones entail flat (or, possibly, downward) curves. In the third row, the weights of the *curvature* factor are presented. The name of this factor follows from the fact that the weights have a more pronounced effect on the middle range maturities of the curve. These concepts are graphically illustrated below.  

```{r, fig.cap = "Yield loadings on the spanned factors", echo=FALSE}
LabSpaFac <- c("Level", "Slope", "Curvature")
N <- length(LabSpaFac)
 
w_pca <- data.frame(t(w[1:N,]))
colnames(w_pca) <- LabSpaFac
w_pca$mat <- mat

# Prepare plots
colors <- c("Level" = "blue", "Slope" = "green", "Curvature" = "red")

g <-  ggplot2::ggplot(data = w_pca, ggplot2::aes(x=  mat)) +  
    ggplot2::geom_line(ggplot2::aes(y = Level, color = "Level"), size = 0.7) + 
    ggplot2::geom_line(ggplot2::aes(y = Slope, color = "Slope"), size = 0.7) +
    ggplot2::geom_line(ggplot2::aes(y = Curvature, color = "Curvature"),  size = 0.7) +
      ggplot2::labs(color = "Legend") + ggplot2::scale_color_manual(values = colors) + ggplot2::theme_classic() +
    ggplot2::theme(axis.title.y= ggplot2::element_blank(), legend.position="top", legend.title=ggplot2::element_blank(), legend.text= ggplot2::element_text(size=8) ) + 
   ggplot2::xlab("Maturity (Years)") + ggplot2::geom_hline(yintercept=0)

print(g)
```

To directly obtain the time series of the country-specific spanned factors, the user can simply use the *Spanned_Factors* function as follows:

```{r}
data('CM_Yields')
Economies <- c("China", "Brazil", "Mexico", "Uruguay")
N <- 2
SpaFact <- Spanned_Factors(Yields, Economies, N)
```


## The P-dynamics estimation
As presented in Table \@ref(tab:ModFea), the dynamics of the risk factors under the $P$-measure evolve according to a VAR(1) model, which may be fully unrestricted (as in the case of the JPS-related models) or somewhat restricted (as in the GVAR and JLL frameworks). Below, the usage of each of these dimensions are illustrated.


### VAR {#S:SectionVAR}
Using the *VAR* function of this package requires simply selecting the appropriated set of risk factors for the desired estimated model. For instance, for the models *JPS jointP* and *VAR jointQ*, the estimation of the $P$-dynamics parameters is obtained as  
```{r}
data("CM_Factors")
PdynPara <- VAR(RiskFactors, VARtype= "unconstrained")
```
whereas the estimation of a *JPS* model for China is 
```{r}
FactorsChina <- RiskFactors[1:7,]
PdynPara <- VAR(FactorsChina, VARtype= "unconstrained")
```

In both cases, the outputs generated are the vector of intercepts and both the feedback and variance-covariance matrices.


### GVAR {#S:SectionGVAR}
The estimation of a GVAR model requires defining a measure of interdependence among the countries of the economic system. This information is reported in the transition matrix, the entries of which reflect the degree of interconnection of two entities of this same economic system. In this package, the illustration of the transition matrix is based on the average of the cross-border trade flow weights for the period spanning the years from 2006 to 2019. Note that each row sums up to 1.

```{r}
data("CM_Trade")
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, DataPath = NULL, TradeFlows)
```

Having defined the form of the transition matrix, one can complete the *GVARinputs* variable along the lines discussed in Section \@ref(S:SectionGVARJLLinputs). 
```{r}
data("CM_Factors_GVAR")

GVARinputs <- list()
GVARinputs$Economies <- Economies
GVARinputs$GVARFactors <- FactorsGVAR
GVARinputs$VARXtype <- "unconstrained"
GVARinputs$Wgvar <- W_gvar

N <- 3

GVARpara <- GVAR(GVARinputs, N)
```

A separate routine is provided for computing the foreign-specific factors (also commonly referred as the star variables) used in the estimation of the VARX models.
```{r}
data('CM_Factors')
StaFac <- StarFactors(RiskFactors, Economies, W_gvar)
```

### JLL {#S:SectionJLL}
Calculating the $P$-dynamics parameters in the form proposed by JLL (2015) requires the following inputs: *(i)* the time series of the risk factors in non-orthogonalized form; *(ii)* the number of country-specific spanned factors and *(iii)* the specification of the *JLLinputs* as presented in Section \@ref(S:SectionGVARJLLinputs). See, for instance:

```{r, eval=FALSE}
data("CM_Factors")
N <- 3
JLLinputs <- list()
ModelType <- "JLL original"  
JLLinputs$Economies <- Economies
JLLinputs$DomUnit <- "China"
JLLinputs$WishSigmas <- 1 
JLLinputs$SigmaNonOrtho <- NULL
JLLinputs$JLLModelType <- ModelType
JLLpara <- JLL(RiskFactors, N, JLLinputs)
```

## ATSM estimation
The estimation of ATSM involves defining the model inputs from the log-likelihood function and, subsequently, the parameters used in its optimization process. The structure proposed in this part of the code is in a great extent based on the term structure package by @LeSingleton2018.

### The log-likelihood function 
For the sake of simplicity, one illustrates the construction of a log-likelihood function based on a *VAR jointQ* model. For both the JLL and GVAR-based models, one must further specify the JLLinputs or GVARinputs along the lines described in the Section \@ref(S:SectionGVARJLLinputs).
```{r}
# Inputs to be specified by the user
data("CM_Yields")
data("CM_Factors")
ModelType <- "VAR jointQ"
Yields <- Yields
ZZ <- RiskFactors
Economies <- c("China", "Brazil", "Mexico", "Uruguay") 
mat <-  Maturities(Yields, Economies, UnitYields = "Month")
DataFrequency <- "Monthly"
GlobalVar <- c("GBC", "CPI_OECD") # Global variables
DomVar <- c("Eco_Act", "Inflation") # Domestic variables 
N <- 3 # Number of country-specific spanned factors

# Generate the "Factor Labels" list (necessary preliminarily step) 
FactorLabels <- LabFac(N, DomVar, GlobalVar, Economies, ModelType) 

# Prepare the inputs for log-Likelihood function
ATSMInputs <- InputsForMLEdensity(ModelType, Yields, ZZ, FactorLabels, mat, Economies, DataFrequency)

# Log-Likelihood function
f <- Functionf(ATSMInputs, Economies, mat, DataFrequency, FactorLabels, ModelType)
```


### Optimization parameters
JPS (2014) requires the estimation of a set of parameters containing 6 elements, namely: the risk-neutral
long-run mean of the short rate (*r0*), the risk-neutral feedback matrix (*K1XQ*), the variance-covariance matrix (*SSZ*) from the VAR processes, the standard deviation of the errors from the portfolios of yields observed with error (*se*), in addition to the intercept (*K0Z*) and the feedback (*K1Z*) matrices of the physical dynamics. Each of these parameters must be cast in an individual list that should contain the *(i)* the starting value of the parameter (if any); *(ii)* the variable label followed by a ':' and a type of constraint; *(iii)* a lower bound (if any), and *(iv)* an upper bound (if any). The variable labels of *r0*, *se*, *K0Z* and *K1Z* should be preceded by the symbol @ as a manner to account for the fact that the solution of these parameters are known in closed-form. For the remaining ones (*K1XQ* and *SSZ*), @LeSingleton2018 provide standardized routines (already incorporated in this package) to set good initial values so that the optimization process runs fast.

The type of constraint to be set for the parameters *r0*, *se*, *K0Z* and *K1Z* is typically *bounded* and are used for bounded matrices. For *K1XQ* and *SSZ*, the *MultiATSM* package currently provides the following options:

- **K1XQ:** *'Jordan'* (or *Jordan MultiCountry'*) for a matrix of Jordan type for single (multi) country specifications. These labels can be extended by *'; stationary'*, if one wishes to impose the largest eigenvalue of the risk-neutral feedback matrix to be strictly smaller than 1;   
- **SSZ:** *'psd'* for a positive semi-definite matrix; *'diag'* for a diagonal matrix; *'BlockDiag'* for a block diagonal matrix (typical from the GVAR-based models) and *'JLLstructure'* for the models containing the restrictions along the lines of JLL (2015). 

One example for the list of parameters specification is

```{r}
K1XQinputs <- list(NULL, "K1XQ: Jordan" , NULL , NULL)
SSZinputs <- list(NULL, "SSZ: psd", NULL, NULL)
r0inputs <- list(NULL, "@r0: bounded", NULL, NULL)
seinputs <- list(NULL, "@se: bounded", 1e-6, NULL)
K0Zinputs<- list(NULL, "@K0Z: bounded", NULL, NULL)
K1Zinputs<- list(NULL, "@K1Z: bounded", NULL, NULL)
```


To complete the optimization setting, the user needs to specify the level of convergence tolerance (usually, 1e-4 is a reasonable value) and the vector *OptRun*. The first element of this vector must be filled by the word *'iter off'*, if the user wishes to switch off the printouts of the numerical optimization routines, or simply by *'iter'* otherwise. The second element of the vector concerns the algorithm that is used in the optimization: the option *'fminunc only'* only uses fminunc, whereas *'fminsearch only'* only applies fminsearch. If *OptRun* is formed exclusively by *'iter off'*, then both fminunc and fminsearch are employed in the optimization.    


# Examples of full implementation of ATSMs {#S:SectionReplication}
This section presents some examples on how to use the *MultiATSM* package to fully implement ATSMs. Further on, the functions of this package are used to replicate some of the results presented in the original papers of JPS (2014) and CM (2021).  


## General template
```{r, eval=FALSE}
###############################################################################################################
#################################### USER INPUTS ##############################################################
###############################################################################################################
# A) Load database data
data("CM_Factors")
data('CM_Factors_GVAR')
data('CM_Trade')
data('CM_Yields')

# B) Decide on general model inputs
ModelType <- "JLL original" # available options are "JPS", "JPS jointP", "GVAR sepQ", "VAR jointQ", "GVAR jointQ", "JLL original", "JLL NoDomUnit", "JLL jointSigma".

StationarityUnderQ <- 1 # Wish to impose stationary condition for the eigenvalues of each country: YES: 1, NO:0

Economies <- c("China", "Brazil", "Mexico", "Uruguay") # Names of the economies from the economic system
GlobalVar <- c("GBC", "CPI_OECD") # Global Variables
DomVar <- c("Eco_Act", "Inflation") # Country-specific variables

N <- 3 # Number of spanned factors per country

OutputLabel <- "Test" # label of the model for saving the file
DataFrequency <- "Monthly" # Frequency of the data
UnitMatYields <- "Month" # time-unit in which yields are expressed. Available options are "Month" or "Year"

# B.1) Decide on specific model inputs
#################################### GVAR-based models #########################################################
if (ModelType == 'GVAR sepQ' || ModelType == 'GVAR jointQ'){
  # Transition matrix inputs:
  t_First <- "2006" # First year of the sample
  t_Last <-  "2019" # Last year of the sample
  W_type <- 'Sample Mean' # Method to compute the transition matrix
  # VARX input:
  VARXtype <- "unconstrained"
}

#################################### JLL-based models #########################################################
if (ModelType =="JLL original" || ModelType == "JLL NoDomUnit" || ModelType == "JLL jointSigma"){
  DominantCountry <- "China" # (i) one of the countries of the system or (ii) "None"
}
################################################################################################################

# C) Decide on Settings for numerical outputs
Horiz <- 50
DesiredGraphs <- c("IRF", "FEVD") #c("Fit", "IRF", "FEVD", "GIRF", "GFEVD")
WishGraphRiskFac <- 0
WishGraphYields <- 1
WishOrthoJLLgraphs <- 0

# D) Bootstrao settings
WishBootstrap <- 1 #  YES: 1; No = 0.
Bootlist <- list()
Bootlist$methodBS <- 'bs' # (i) 'bs' ; (ii) 'wild'; (iii) 'block'
Bootlist$BlockLength <- 4 # necessary input if one chooses the block bootstrap
Bootlist$ndraws <-  30
Bootlist$pctg   <-  95 # confidence level

# E) Out-of-sample forecast
WishForecast <- 1 #  YES: 1; No = 0.
ForecastList <- list()
ForecastList$ForHoriz <- 12 # forecast horizon
ForecastList$t0Sample <- 1 # initial sample date
ForecastList$t0Forecast <- 145 # last sample date for the first forecast

###############################################################################################################
############################### NO NEED TO MAKE CHANGES FROM HERE #############################################
###############################################################################################################

# 2) Minor preliminary work
C <- length(Economies)
FactorLabels <- LabFac(N, DomVar,GlobalVar, Economies, ModelType) # Generate the set of labels
mat <- Maturities(Yields, Economies, UnitYields = UnitMatYields) # Vector of maturities expressed in years
# (alternatively, it can be completed manually)
ZZ <- RiskFactors
# 2.1) Build the GVARinputs
if (ModelType == 'GVAR sepQ' || ModelType == 'GVAR jointQ'){
  GVARinputs <- list()
  GVARinputs$Economies <- Economies
  GVARinputs$GVARFactors <- FactorsGVAR
  GVARinputs$VARXtype <- VARXtype
  GVARinputs$Wgvar <- Transition_Matrix(t_First, t_Last, Economies, W_type, DataPath = NULL, Data = TradeFlows)
} else { GVARinputs <- NULL }

# 2.2) Build the JLLinputs
if (ModelType =="JLL original" || ModelType == "JLL NoDomUnit" || ModelType == "JLL jointSigma"){
  JLLinputs <- list()
  JLLinputs$Economies <- Economies
  JLLinputs$DomUnit <- DominantCountry
  JLLinputs$WishSigmas <- 1 # Sigma matrix is estimated within the "InputsForMLEdensity" function
  JLLinputs$SigmaNonOrtho <- NULL
  JLLinputs$JLLModelType <- ModelType
}else{
  JLLinputs <- NULL
}

# 3) Prepare the inputs of the likelihood function
ModelParaList <- list()
for (i in 1:C){
  if (( ModelType == "GVAR jointQ" || ModelType == "VAR jointQ" || ModelType == "JLL original"
        || ModelType == "JLL NoDomUnit" || ModelType == "JLL jointSigma" )   & i >1 ){break} # the break avoids the models with joint estimation under the Q to be estimated C times

  # 3.1) Compute the inputs that go directly into the log-likelihood function
  ATSMInputs <- InputsForMLEdensity(ModelType, Yields, ZZ, FactorLabels, mat, Economies, DataFrequency, JLLinputs,   GVARinputs)

  # 3.2) Initial guesses for Variables that will be concentrared out of from the log-likelihood function
  K1XQ <- ATSMInputs$K1XQ
  if (ModelType == "JLL original" || ModelType == "JLL NoDomUnit" ){ SSZ <- NULL} else{SSZ <- ATSMInputs$SSZ}

  # 4) Build the objective function
  f <- Functionf(ATSMInputs, Economies, mat, DataFrequency, FactorLabels, ModelType)

  # 5) Set the optimization settings
  VarLab <- ParaLabels(ModelType, StationarityUnderQ)

  varargin <- list()
  varargin$K1XQ <-list(K1XQ, VarLab[[ModelType]][["K1XQ"]] , NULL , NULL)
  varargin$SSZ <- list(SSZ, VarLab[[ModelType]][["SSZ"]], NULL, NULL)
  varargin$r0 <- list(NULL, VarLab[[ModelType]][["r0"]], NULL, NULL)
  varargin$se <- list(NULL, VarLab[[ModelType]][["se"]], 1e-6, NULL)
  varargin$K0Z <- list(NULL, VarLab[[ModelType]][["K0Z"]], NULL, NULL)
  varargin$K1Z <- list(NULL, VarLab[[ModelType]][["K1Z"]], NULL, NULL)
  varargin$OptRun <- c("iter off")

  LabelVar<- c('Value', 'Label', 'LB', 'UB') # Elements of each parameter
  for (d in 1:(length(varargin)-1)){ names(varargin[[d]]) <-  LabelVar}

  tol <- 1e-4

  # 6) Optimization of the model
  if (ModelType == 'JPS' || ModelType == 'JPS jointP' || ModelType == "GVAR sepQ"){
    ModelParaList[[ModelType]][[Economies[i]]] <- Optimization(f, tol, varargin, FactorLabels, Economies, ModelType, JLLinputs, GVARinputs)$Summary
  }else{ ModelParaList[[ModelType]] <- Optimization(f, tol, varargin, FactorLabels, Economies, ModelType, JLLinputs, GVARinputs)$Summary}

}

# 7) Numerical and graphical outputs
InputsForOutputs <- InputsForOutputs(ModelType, Horiz, DesiredGraphs, OutputLabel, StationarityUnderQ,                  WishGraphYields, WishGraphRiskFac, WishOrthoJLLgraphs, WishBootstrap, Bootlist, WishForecast, ForecastList)
# A) Fit, IRF, FEVD, GIRF and GFEVD
NumericalOutputs <- NumOutputs(ModelType, ModelParaList, InputsForOutputs, FactorLabels, Economies)

# B) Bootstrap
Bootstrap <- Bootstrap(ModelType, ModelParaList, NumericalOutputs, mat, Economies, InputsForOutputs, FactorLabels, DataFrequency, varargin, JLLinputs, GVARinputs)

# C) Out-of-sample forecasting
Forecasts <- ForecastYields(ModelType, ModelParaList, InputsForOutputs, FactorLabels, Economies, DataFrequency, JLLinputs, GVARinputs)

```



## Paper replications
### Joslin, Priebisch and Singleton (2014)
To the best of my knowledge, the original data-set used in @JoslinPriebschSingleton2014 is not publicly available. Nevertheless, there exists an alternative source of data that closely matches JPS (2014) in terms of the construction method and the sample period. This database is constructed by @BauerRudebusch2017 (BR, 2017) and is available at [Bauer's website](https://www.michaeldbauer.com/publication/spanning-puzzle/). In that paper, BR (2017) investigate whether macro-finance term structure models better suit the unspanned macro risk framework by JPS (2014) or the earlier traditional spanned settings as the one by @AngPiazzesi2003JME. Accordingly, BR (2017) intend to replicate some of the empirical results reported in JPS (2014). The R-code used by BR (2017) is also available at [Bauer's website](https://www.michaeldbauer.com/publication/spanning-puzzle/).

Supported by BR (2017)'s data-set, the code below uses the *MultiATSM* package to estimate the key model parameters from the ATSM along the lines of JPS (2014).

```{r, eval=FALSE} 
###############################################################################################################
#################################### USER INPUTS ##############################################################
###############################################################################################################
# A) Load database data
data("BR_jps_gro_R3")

# B) Decide on general model inputs
ModelType <- "JPS"

StationarityUnderQ <- 0 # Wish to impose stationary condition for the eigenvalues of each country: YES: 1, NO:0

Economies <- c("U.S.") # Names of the economies from the economic system
GlobalVar <- c()
DomVar <- c("GRO", "INF") # Country-specific variables

N <- 3 # Number of spanned factors per country

DataFrequency <- "Monthly" # Frequency of the data
UnitMatYields <- "Month" # time-unit in which yields are expressed. Available options are "Month" or "Year"

mat <- c(0.25, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10) # maturities expressed in years

###############################################################################################################
############################### NO NEED TO MAKE CHANGES FROM HERE #############################################
###############################################################################################################
# 2) Minor preliminary work
C <- length(Economies)
M <- length(DomVar)
FactorLabels <- LabFac(N, DomVar, GlobalVar, Economies, ModelType) # Generate the set of labels
Yields <- t(BR_jps_out$Y)
W <- BR_jps_out$W[1:N,] # Use the Wpca matrix from BR (2017)
SpaFac <- W%*%Yields
rownames(SpaFac) <- FactorLabels$Tables$U.S.[-(1:M)]
ZZ <- rbind(t(BR_jps_out$M.o), SpaFac) # Complete set of risk factors


# 3) Prepare the inputs of the likelihood function
for (i in 1:C){
# 3.1) Compute the inputs that go directly into the log-likelihood function
ATSMInputs <- InputsForMLEdensity(ModelType, Yields, ZZ, FactorLabels, mat, Economies, DataFrequency)
# 3.2) Initial guesses for Variables that will be concentrared out of from the log-likelihood function
K1XQ <- ATSMInputs$K1XQ
SSZ <- ATSMInputs$SSZ
# 3.3) Adjust the inputs which are funtion of the W matrix
ATSMInputs$Wpca <- W
ATSMInputs$We <- t(pracma::null(W))
ATSMInputs$WpcaFull <- rbind(ATSMInputs$Wpca, ATSMInputs$We)
ATSMInputs$PP <- SpaFac

# 4) Build the objective function
f <- Functionf(ATSMInputs, Economies, mat, DataFrequency, FactorLabels, ModelType)

# 5) Set the optimization settings
VarLab <- ParaLabels(ModelType, StationarityUnderQ)

varargin <- list()
varargin$K1XQ <-list(K1XQ, VarLab[[ModelType]][["K1XQ"]] , NULL , NULL)
varargin$SSZ <- list(SSZ, VarLab[[ModelType]][["SSZ"]], NULL, NULL)
varargin$r0 <- list(NULL, VarLab[[ModelType]][["r0"]], NULL, NULL)
varargin$se <- list(NULL, VarLab[[ModelType]][["se"]], 1e-6, NULL)
varargin$K0Z <- list(NULL, VarLab[[ModelType]][["K0Z"]], NULL, NULL)
varargin$K1Z <- list(NULL, VarLab[[ModelType]][["K1Z"]], NULL, NULL)
varargin$OptRun <-  c("iter off")

LabelVar<- c('Value', 'Label', 'LB', 'UB') # Elements of each parameter
for (d in 1:(length(varargin)-1)){ names(varargin[[d]]) <-  LabelVar}

tol <- 1e-4

# 6) Optimization of the model
ModelPara <- Optimization(f, tol, varargin, FactorLabels, Economies, ModelType)$Summary
}
```



The tables below compare the ATSM parameter estimates generated from BR (2017) and the *MultiATSM*. Overall, one can note that the differences in the estimates are economically modest.


```{r, echo= FALSE}
options(scipen = 100) # eliminate the scientific notation
data("BR_jps_gro_R3")
data("JPSrep")

RowsQ <- c("$r0$", "$\\lambda_1$", "$\\lambda_2$", "$\\lambda_3$" )
TableQ <- data.frame(matrix(NA, ncol = 0, nrow =length(RowsQ)))
row.names(TableQ) <- RowsQ

PackageQ<- c(ModelPara$ests$r0, diag(ModelPara$ests$K1XQ))
BRq <- c(BR_jps_out$est.llk$rho0.cP, diag(BR_jps_out$est.llk$KQ.XX))
TableQ$MultiATSM <- PackageQ
TableQ$'BR (2017)' <- BRq

TableQ <- round(TableQ, digits = 5)


kableExtra::kbl(TableQ, align = "c", caption = "$Q$-dynamics parameters") %>%
  kableExtra::kable_classic("striped", full_width = F)  %>%
  kableExtra::row_spec(0, font_size = 14) %>%
  kableExtra::footnote(general = " $\\lambda$'s are the eigenvalues from the risk-neutral feedback matrix and $r0$ is the long-run mean of the short rate under Q.")
```


```{r, echo= FALSE}
data("BR_jps_gro_R3")
data("JPSrep")

RowsP <- c("PC1", "PC2", "PC3", "GRO", "INF")
ColP <- c(" ", RowsP)

# 1) K0Z and K1Z
# Bauer and Rudebusch coefficients
TablePbr <- data.frame(matrix(NA, ncol = length(ColP), nrow =length(RowsP)))
row.names(TablePbr) <- RowsP
colnames(TablePbr) <- ColP

TablePbr[[ColP[1]]] <- BR_jps_out$est.llk$KP.0Z
for(j in 1:length(RowsP) ){TablePbr[[RowsP[j]]] <- BR_jps_out$est.llk$KP.ZZ[,j]}

TablePbr <- round(TablePbr, digits = 5)

# MultiATSM coefficients
TablePMultiATSM <- data.frame(matrix(NA, ncol = length(ColP), nrow =length(RowsP)))
row.names(TablePMultiATSM) <- RowsP
colnames(TablePMultiATSM) <- ColP

IdxVar <- c(3:5, 1:2) # indexes to flip the order of the spanned and unspanned factors
TablePMultiATSM[[ColP[1]]] <- ModelPara$ests$K0Z[IdxVar]
ModelPara$ests$K1Z <- ModelPara$ests$K1Z[IdxVar,IdxVar]
for(j in 1:length(RowsP) ){ TablePMultiATSM[[RowsP[j]]] <- ModelPara$ests$K1Z[,j]}


TablePMultiATSM <- round(TablePMultiATSM, digits = 5)
TableP <- rbind(TablePbr,TablePMultiATSM)
row.names(TableP) <- c(RowsP,paste(RowsP," ",sep="")) # trick to avoid rows to have the same name

kableExtra::kbl(TableP, align = "c", caption = "$P$-dynamics parameters") %>%
  kableExtra::kable_classic("striped", full_width = F)  %>%
  kableExtra::row_spec(0, font_size = 14) %>%
  kableExtra::add_header_above(c(" "= 1, "K0Z"=1, "K1Z" = 5), bold = T) %>%
  kableExtra::pack_rows("BR (2017)", 1, 5) %>%
  kableExtra::pack_rows("MultiATSM", 6, 10) %>%
  kableExtra::footnote(general = " $K0Z$ is the intercept and $K1Z$ is feedback matrix from the $P$-dynamics.")
```

```{r, echo= FALSE}
data("BR_jps_gro_R3")
data("JPSrep")

se <- data.frame(BR_jps_out$est.llk$sigma.e, ModelPara$ests$se )
rownames(se) <- "se"
colnames(se) <- c("MultiATSM","BR (2017)")
se <- round(se, digits = 7)

kableExtra::kbl(se, align = "c", caption ="Portfolio of yields with errors") %>%
  kableExtra::kable_classic("striped", full_width = F)  %>%
  kableExtra::row_spec(0, font_size = 14) %>%
  kableExtra::footnote(general = " $se$ is the standard deviation of the portfolio of yields observed with errors.")
```


It is relevant to highlight that the script above makes use of 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 this 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 form of the log-likelihood function (llk): while in the BR (2017) code the llk 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 llk up to the Jacobian term).



### Candelon and Moura (2021)
```{r, eval=FALSE}
###############################################################################################################
#################################### USER INPUTS ##############################################################
###############################################################################################################
# A) Load database data
data("CM_Factors")
data('CM_Factors_GVAR')
data('CM_Trade')
data('CM_Yields')

# B) Decide on general model inputs
ModelType <- "GVAR jointQ" # Set "GVAR jointQ" for the GVAR-ATSM and "JLL original" for the JLL-ATSM

StationarityUnderQ <- 0

Economies <- c("China", "Brazil", "Mexico", "Uruguay") 
GlobalVar <- c("GBC", "CPI_OECD") # Global Variables
DomVar <- c("Eco_Act", "Inflation") # Country-specific variables

N <- 3 # Number of spanned factors per country

OutputLabel <- "CM_2021" # label of the model for saving the file
DataFrequency <- "Monthly" # Frequency of the data
UnitMatYields <- "Month" # time-unit in which yields are expressed. Available options are "Month" or "Year"

# B.1) Decide on specific model inputs
#################################### GVAR-based models #########################################################
if (ModelType == 'GVAR jointQ'){
# Transition matrix inputs:
t_First <- "2006" # First year of the sample
t_Last <-  "2019" # Last year of the sample
W_type <- 'Sample Mean' # Method to compute the transition matrix
# VARX input:
VARXtype <- "unconstrained"
}

#################################### JLL-based models #########################################################
if (ModelType =="JLL original"){
    DominantCountry <- "China"
}
################################################################################################################

# C) Decide on Settings for numerical outputs
Horiz <- 25
DesiredGraphs <- c("Fit", "GIRF", "GFEVD")
WishGraphRiskFac <- 0
WishGraphYields <- 1
WishOrthoJLLgraphs <- 1

# D) Bootstrao settings
WishBootstrap <- 0 #  If bootstrap outputs are desired set this variable equal to 1 (it may take several hours/days).
Bootlist <- list()
Bootlist$methodBS <- 'bs'
Bootlist$BlockLength <- c()
Bootlist$ndraws <-  1000
Bootlist$pctg   <-  95

# E) Out-of-sample forecast
WishForecast <- 0 #  If out-of-sample forecasts outputs are desired set this variable equal to 1
ForecastList <- list()
ForecastList$ForHoriz <- 12
ForecastList$t0Sample <- 1
ForecastList$t0Forecast <- 90

###############################################################################################################
############################### NO NEED TO MAKE CHANGES FROM HERE #############################################
###############################################################################################################

# 2) Minor preliminary work
C <- length(Economies)
FactorLabels <- LabFac(N, DomVar,GlobalVar, Economies, ModelType) # Generate the set of labels
mat <- Maturities(Yields, Economies, UnitYields = UnitMatYields) # Vector of maturities expressed in years

ZZ <- RiskFactors
# 2.1) Build the GVARinputs
if (ModelType == "GVAR jointQ"){
  GVARinputs <- list()
  GVARinputs$Economies <- Economies
  GVARinputs$GVARFactors <- FactorsGVAR
  GVARinputs$VARXtype <- VARXtype
  GVARinputs$Wgvar <- Transition_Matrix(t_First, t_Last, Economies, W_type, DataPath = NULL, Data = TradeFlows)
  } else { GVARinputs <- NULL }

# 2.2) Build the JLLinputs
if (ModelType =="JLL original"){
  JLLinputs <- list()
  JLLinputs$Economies <- Economies
  JLLinputs$DomUnit <- DominantCountry
  JLLinputs$WishSigmas <- 1
  JLLinputs$SigmaNonOrtho <- NULL
  JLLinputs$JLLModelType <- ModelType
}else{  JLLinputs <- NULL }

# 3) Prepare the inputs of the likelihood function
ModelParaList <- list()

# 3.1) Compute the inputs that go directly into the log-likelihood function
ATSMInputs <- InputsForMLEdensity(ModelType, Yields, ZZ, FactorLabels, mat, Economies, DataFrequency, JLLinputs,
GVARinputs)

  # 3.2) Initial guesses for Variables that will be concentrared out of from the log-likelihood function
K1XQ <- ATSMInputs$K1XQ
if (ModelType == "JLL original"){ SSZ <- NULL} else{SSZ <- ATSMInputs$SSZ}

# 4) Build the objective function
f <- Functionf(ATSMInputs, Economies, mat, DataFrequency, FactorLabels, ModelType)

# 5) Set the optimization settings
VarLab <- ParaLabels(ModelType, StationarityUnderQ)

varargin <- list()
varargin$K1XQ <-list(K1XQ, VarLab[[ModelType]][["K1XQ"]], NULL , NULL)
varargin$SSZ <- list(SSZ, VarLab[[ModelType]][["SSZ"]], NULL, NULL)
varargin$r0 <- list(NULL, VarLab[[ModelType]][["r0"]], NULL, NULL)
varargin$se <- list(NULL, VarLab[[ModelType]][["se"]], 1e-6, NULL)
varargin$K0Z <- list(NULL, VarLab[[ModelType]][["K0Z"]], NULL, NULL)
varargin$K1Z <- list(NULL, VarLab[[ModelType]][["K1Z"]], NULL, NULL)
varargin$OptRun <-  c("iter off")

LabelVar<- c('Value', 'Label', 'LB', 'UB')
for (d in 1:(length(varargin)-1)){ names(varargin[[d]]) <-  LabelVar}

tol <- 1e-4

# 6) Optimization of the model
ModelParaList[[ModelType]] <- Optimization(f, tol, varargin, FactorLabels, Economies, ModelType, JLLinputs, GVARinputs)$Summary


# 7) Numerical and graphical outputs
InputsForOutputs <- InputsForOutputs(ModelType, Horiz, DesiredGraphs, OutputLabel, StationarityUnderQ, WishGraphYields,  WishGraphRiskFac, WishOrthoJLLgraphs, WishBootstrap, Bootlist, WishForecast, ForecastList)

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

# B) Bootstrap
Bootstrap <- Bootstrap(ModelType, ModelParaList, NumericalOutputs, mat, Economies, InputsForOutputs, FactorLabels,
DataFrequency, varargin, JLLinputs, GVARinputs)

# C) Out-of-sample forecasting
Forecasts <- ForecastYields(ModelType, ModelParaList, InputsForOutputs, FactorLabels, Economies, DataFrequency, JLLinputs, GVARinputs)
```



# References

