---
title: "Replication file for: Anchoring the yield curve: term structure forecasting using survey data"
subtitle: "DNS model in state space form"
date: \today
author:
- Carlo Altavilla^[European Central Bank]
- Raffaella Giacomini^[University College London]
- Giuseppe Ragusa^[Luiss University]
output:
  pdf_document:
    keep_tex: true
    includes:
      in_header: header.tex
    toc: false
fontsize: 12pt
geometry: margin=1in
---

\listoffigures
\listoftables
\clearpage

```{r load_packages, include=TRUE}
if (!require("pacman")) {
  install.packages("pacman", repos = "https://cran.rstudio.com")
  require(pacman)
  }

pkg_from_cran = c(
  "tikzDevice",
  "Hmisc",
  "RColorBrewer",
  "RcppRoll",
  "dplyr",
  "forecast",
  "ggplot2",
  "ggthemes",
  "lubridate",
  "purrr",
  "reshape2",
  "sandwich",
  "xts"
  )

p_load(char = pkg_from_cran)

## Patch latex form Hmisc to work with markdown
if (!is.na(match("latex.R", list.files("."))))
  source("latex.R")

```


```{r knitr_setup, include=FALSE, cache = FALSE}
library(knitr)
a4width  <- 8.3
a4height <- 11.7

opts_chunk$set(
  cache = FALSE,
  dev = 'tikz',
  size = 'footnotesize',
  fig.path = 'figures/',
  cache.path = 'cache-anchored/',
  fig.width = a4width*.65,
  fig.height = a4height*.3,
  fig.align = 'center',
  error = FALSE,
  tidy = TRUE,
  message = FALSE,
  echo = FALSE
  )

options(width = 55, formatR.arrow = TRUE)
options(stringsAsFactors = FALSE) ## PITA


```





```{r}
cap_0 <- "DNS factors and empirical counterparts."
note_0 <- "The first factor $\\beta_{1t}$ controls the yield curve level, as it can be verified that $\\lim_{\\tau\\to\\infty}y_{t}(\\tau)=\\beta_{1t}$. The second factor $\\beta_{2t}$ is related to the yield curve slope, defined as the difference between the 10-year and three-month yields. The third factor $\\beta_{3t}$ governs the curvature of the yield curve, defined as twice the two-year yield minus the sum of the 10-year and three-month yields."
```


```{r external-data, echo=FALSE, message=FALSE, warning=FALSE}
load("../results/dns_kalman.Rda")
source("../R/01-helper.R")
source("../R/02-load_data.R")
source("functions.R")
```

```{r survey_rel_msfe, results="asis"}
off <- which(index(yield_data) == "Jan 2000")
bc_error <- list()
mat_ <- c(1, 2, 4, 8, 12, 17)
excl <- -c(1:(off - 1))

for (mm in c(3, 6, 12, 24, 60, 120)) {
  mat <- paste0("m", mm)
  mu  <- bc_set[[mat]][[2]]
  idx <- match(mm, maturities)
  for (hh in c(3, 6, 9, 12)) {
    bc_error[[mat]][[paste0("h=", hh)]] <-
      yield_data[excl, idx] - mu[[hh / 3]][excl]
  }
}
msfe_bc <- matrix(0, 6, 4)
for (mm in 1:6) {
  msfe_bc[mm, ] <- sapply(bc_error[[mm]], MSFE)
}

msfe_rw <- map(dns_model, ~ MSFE(.x$error_rw))[[1]][seq(3, 12, 3), c(1, 2, 4, 8, 12, 17)]

tbl <- t(t(msfe_bc)/msfe_rw)
colnames(tbl) <- paste0("$h = ", seq(3,12,3), "$")
rownames(tbl) <- paste0("$\\tau = ", c(3, 6, 12, 25, 60, 120), "$")

latex.def(
  format(tbl, digits = 2),
  dcolumn = TRUE,
  file = "",
  title = "",
  rowlabel = "Maturity",
  caption = "MSFE of Blue Chip survey forecasts (relative to the Random Walk forecast)",
  caption.loc = "bottom",
  nomargin = FALSE,
  numeric.dollar = FALSE,
  comment = FALSE
  )

```

```{r survey_rel_msfe_save, results="hide"}
w <- latex.def(
  format(t(tbl), digits = 2),
  dcolumn = TRUE,
  file = "tables/survey_rel_msfe.tex",
  title = "",
  rowlabel = "Maturity",
  rownames = rownames(t(tbl)),
  caption = "MSFE of Blue Chip survey forecasts (relative to the Random Walk forecast)",
  caption.loc = "bottom",
  nomargin = FALSE,
  table.env = FALSE,
  numeric.dollar = FALSE,
  comment = FALSE
  )
```


```{r survey_rel_dns_msfe, results="asis"}


msfe_dns <- map(dns_model, ~ MSFE(.x$error))[[1]][seq(3, 12, 3), c(1, 2, 4, 8, 12, 17)]

tbl <- t(t(msfe_bc)/msfe_dns)
colnames(tbl) <- paste0("$h = ", seq(3,12,3), "$")
rownames(tbl) <- paste0("$\\tau = ", c(3, 6, 12, 24, 60, 120), "$")

latex.def(
  format(tbl, digits = 2),
  dcolumn = TRUE,
  file = "",
  title = "",
  rowlabel = "Maturity",
  caption = "MSFE of Blue Chip survey forecasts (relative to the Yields-only DNS  forecasts)",
  caption.loc = "bottom",
  nomargin = FALSE,
  numeric.dollar = FALSE,
  comment = FALSE
  )

```

```{r survey_rel_dns_msfe_save, results="hide"}
w <- latex.def(
  format(tbl, digits = 2),
  dcolumn = TRUE,
  file = "tables/survey_rel_dns_msfe.tex",
  title = "",
  rowlabel = "Maturity",
  #rownames = rownames(tbl),
  caption = "MSFE of Blue Chip survey forecasts (relative to the Yields-only DNS forecasts).",
  caption.loc = "bottom",
  nomargin = FALSE,
  table.env = FALSE,
  numeric.dollar = FALSE,
  comment = FALSE,

  )
```





```{r factors, fig.cap=paste0(cap_0, note_0)}
emp_level <- (yield_data$'3' + yield_data$'24' + yield_data$'120')/3
emp_slope <- (yield_data$'3' - yield_data$'120')
emp_curva <- 2*yield_data$'24' - yield_data$'120' - yield_data$'3'

par(mfcol = c(3, 1), mar = c(2, 3, 2, 1), mgp = c(2, .7, 0), tck = -.01)

mm <- MARSS::MARSSkfas(dns_model[[1]]$est_model[[145]], return.kfas.model = TRUE)

BX <- cbind(mm$xtT[1,], emp_level[1:324], mm$xtT[2,], emp_slope[1:324], mm$xtT[3,], emp_curva[1:324])
BX <- ts(as.matrix(BX), start = c(1985,1), end = c(2011, 12), frequency = 12)

plot(BX[,1], ylim = c(0,12), ylab = "Percent", lwd = 2, xlab = "")
lines(BX[,2], lty = 3, lwd = 2)
grid(lwd = 1.2)
legend(x = 'topright', legend = c('Level factor', '(y(3)+y(24)+y(120))/3'), lty =
         c(1, 3), bty = 'n')


plot(BX[,3], lwd = 2, ylab = 'Percent', ylim = c(-8,4), xlab = "")
lines(BX[,4], lty = 3, lwd = 2)
grid(lwd = 1.2)
legend(
  x = 'topright',
  legend = c('Slope factor', 'y(3)-y(120)'),
  lty = c(1, 3),
  bty = 'n'
  )


plot(0.37*BX[,5], lwd = 2, ylab = 'Percent', xlab = "", ylim = c(-4, 2))
lines(BX[, 6], lty = 3, lwd = 2)
grid(lwd = 1.2)
legend(x = 'topright', legend = c('0.37*Curvature factor', '2 * y(24)-y(120)-y(3)'),
       lty = c(1, 3), bty = 'n')

```



```{r insmplfit, results='asis'}
mm <- MARSS::MARSSkfas(dns_model[[1]]$est_model[[145]], return.kfas.model = TRUE)
yhat <- as.matrix(fitted(mm$kfas.model))
y <- as.matrix(yield_data[1:324,])
resid_hat <- y - yhat
stat <- apply(resid_hat, 2, function(u) c(mean(u), sd(u), min(u), max(u), MAE(u), RMSFE(u), Acf(u, lag.max = 30, plot = FALSE)$acf[c(2, 13, 31)]))
stat <- t(stat)
stat <- data.frame(stat)
colnames(stat) <- c("Mean", "SD", "Min", "Max", "MAE", "RMSE", "$\\hat{\\rho}(1)$","$\\hat{\\rho}(12)$", "$\\hat{\\rho}(30)$")
# colnames(stat) <- c('Mean', 'SD', 'Min', 'Max', 'MAE', 'RMSE', '$\\rho(1)$','$\\rho(12)$', '$\\rho(30)$')
stat <- cbind(as.numeric(rownames(stat)), stat)
colnames(stat)[1] <- 'Maturity (months)'

latex.def(format(stat, digits = 2), dcolumn = TRUE,  file = "", title = "", rowname = NULL,
      numeric.dollar = FALSE,  size = "footnotesize", caption = "In-sample fit statistic of the DNS model in state-space form. The descriptive statistics refer to the residuals of the equation for yields at the corresponding maturities. The last three columns present residual sample autocorrelations at lag 1, 12, and 30, respectively.", caption.loc = "bottom", nomargin = FALSE, comment = FALSE)
```

```{r insmpfit_save, results='hide'}
w <-
  latex.def(
  format(stat, digits = 2),
  rowname = NULL,
  dcolumn = TRUE,
  file = "tables/insmpfit.tex",
  title = "",
  numeric.dollar = FALSE,
  size = "footnotesize",
  table.env = FALSE,
  nomargin = FALSE,
  comment = FALSE
  )

```

```{r levslpcrv_info, results='hide'}

mean_curv_diff <- mean(0.37*mm$xtT[3,] - emp_curva[1:324])
sd_curv_diff <- sd(0.37*mm$xtT[3,] - emp_curva[1:324])

mean_slope_diff <- mean(mm$xtT[2,] - emp_slope[1:324])
sd_slope_diff <- sd(mm$xtT[2,] - emp_slope[1:324])
cor_slope_diff <- cor(mm$xtT[2,], emp_slope[1:324])


mean_level_diff <-  mean(mm$xtT[1,] - emp_level[1:324])
sd_level_diff <- sd(mm$xtT[1,] - emp_level[1:324])

mean_level_diff_upDec01 <-  mean(mm$xtT[1,1:204] - emp_level[1:204])
sd_level_diff_upDec01 <- sd(mm$xtT[1, 1:204] - emp_level[1:204])

mean_level_diff_fromJan02 <-  mean(mm$xtT[1,-c(1:204)] - emp_level[205:324])
sd_level_diff_fromJan02 <- sd(mm$xtT[1, -c(1:204)] - emp_level[205:324])


cor_level_diff_upDec01 <- cor(mm$xtT[1,1:204], emp_level[1:204])
cor_level_diff_fromJan02 <- cor(mm$xtT[1, -c(1:204)], emp_level[205:324])
```

The curvature factor closely matches the dynamics of its empirical counterpart: the difference between the two series has a mean of `r mean_curv_diff*100` basis points and a standard deviation of `r sd_curv_diff*100` basis points. Also, the slope factor matches very closely the empirical proxy for the slope with a correlation of `r cor_slope_diff`. The level factor shows instead a marked departure from the empirical counterpart. Importantly, this departure is most noticeable in the period 2000-2011. In particular, from January 1985 to December 2001, the correlation between $\\hat{\beta}_{1,t}$ and $(y_{t}(3)+y_{t}(24)+y_{t}(120))/3$ is `r cor_level_diff_upDec01`, from January 2002 to December 2011 the correlation drops to `r cor_level_diff_fromJan02`. The mean and standard deviations of the difference increase from `r mean_level_diff_upDec01*100` basis points to `r mean_level_diff_fromJan02*100` basis points and from `r sd_level_diff_upDec01*100` basis points to `r sd_level_diff_fromJan02*100` basis points, respectively.



```{r calculate_mse, echo=FALSE}
## dns_model is a list each is a model
maturities <- c(seq(3,24,3), 30, seq(36, 120, 12))
cnames <- paste0('$\\tau=', maturities, "$")
rnames <- paste0('$h=', 1:12, "$")

msfe_ss <- map(dns_model, ~MSFE(.x$error, colnames = cnames, rownames = rnames))
msfe_rw <- map(dns_model, ~MSFE(.x$error_rw, colnames = cnames, rownames = rnames))
msfe_tl <- map(dns_model, ~MSFE(.x$error_tilt, colnames = cnames, rownames = rnames))
msfe_cf <- map(dns_model, ~MSFE(.x$error_cond, colnames = cnames, rownames = rnames))

rel_msfe <- list(
  ss_ss = map2(msfe_ss, msfe_rw, ~ .x / .y),
  ss_tl = map2(msfe_tl, msfe_rw, ~ .x / .y),
  ss_cf = map2(msfe_cf, msfe_rw, ~ .x / .y)
  )

msfe_ss <- map(dns_model, ~MSFE(.x$error[1:108, , ],
                                colnames = cnames, rownames = rnames))
msfe_rw <- map(dns_model, ~MSFE(.x$error_rw[1:108, , ],
                                colnames = cnames, rownames = rnames))
msfe_tl <- map(dns_model, ~MSFE(.x$error_tilt[1:108, , ],
                                colnames = cnames, rownames = rnames))

msfe_tl_2 <- map(dns_model, ~MSFE(.x$error_tilt_2[1:108, , ],
                                colnames = cnames, rownames = rnames))


msfe_cf <- map(dns_model, ~MSFE(.x$error_cond[1:108, , ],
                                colnames = cnames, rownames = rnames))

rel_msfe_up2008 <- list(
  ss_ss = map2(msfe_ss, msfe_rw, ~ .x / .y),
  ss_tl = map2(msfe_tl, msfe_rw, ~ .x / .y),
  ss_tl_2 = map2(msfe_tl_2, msfe_rw, ~ .x / .y),
  ss_cf = map2(msfe_cf, msfe_rw, ~ .x / .y)
  )


msfe_ss <- map(dns_model, ~MSFE(.x$error[49:142, , ],
                                colnames = cnames, rownames = rnames))
msfe_rw <- map(dns_model, ~MSFE(.x$error_rw[49:142, , ],
                                colnames = cnames, rownames = rnames))
msfe_tl <- map(dns_model, ~MSFE(.x$error_tilt[49:142, , ],
                                colnames = cnames, rownames = rnames))

msfe_tl_2 <- map(dns_model, ~MSFE(.x$error_tilt_2[49:142, , ],
                                colnames = cnames, rownames = rnames))

msfe_cf <- map(dns_model, ~MSFE(.x$error_cond[49:142, , ],
                                colnames = cnames, rownames = rnames))

rel_msfe_from2005 <- list(
  ss_ss = map2(msfe_ss, msfe_rw, ~ .x / .y),
  ss_tl = map2(msfe_tl, msfe_rw, ~ .x / .y),
  ss_tl_2 = map2(msfe_tl_2, msfe_rw, ~ .x / .y),
  ss_cf = map2(msfe_cf, msfe_rw, ~ .x / .y)
  )


msfe_ss <- map(dns_model, ~RollMSFE(.x$error, n = 48,
                                    colnames = cnames, rownames = rnames))
msfe_rw <- map(dns_model, ~RollMSFE(.x$error_rw, n = 48,
                                    colnames = cnames, rownames = rnames))
msfe_tl <- map(dns_model, ~RollMSFE(.x$error_tilt, n = 48,
                                    colnames = cnames, rownames = rnames))

msfe_tl_2 <- map(dns_model, ~RollMSFE(.x$error_tilt_2, n = 48,
                                    colnames = cnames, rownames = rnames))


msfe_cf <- map(dns_model, ~RollMSFE(.x$error_cond, n = 48,
                                    colnames = cnames, rownames = rnames))

roll_rel_msfe <- list(
  ss_ss = map2(msfe_ss, msfe_rw, ~ .x / .y),
  ss_tl = map2(msfe_tl, msfe_rw, ~ .x / .y),
  ss_tl_2 = map2(msfe_tl_2, msfe_rw, ~ .x / .y),
  ss_cf = map2(msfe_cf, msfe_rw, ~ .x / .y)
  )


## Create rollin matrix to be plotted using ggplot
ggmat <- c(3, 6, 9, 12, 17)
hh <- c(3, 6, 9, 12)
dat <- seq(ymd("2000-01-2"), ymd("2012-01-2"), by = "months")
roll_mod_info <- list()
for (j in 1:2) {
  XX <- roll_rel_msfe[[1]][[j]][, hh, ggmat]
  XX <- melt(XX)
  XX <- cbind(XX, ifelse(j == 1, "Yields only DNS", "Macro augmented DNS"))
  colnames(XX) <- c("Dates", "h", "tau", "rm", "mm")
  XX$Dates <- dat[-c(1:47)]
  YY <- roll_rel_msfe[[2]][[j]][, hh, ggmat]
  YY <- melt(YY)
  YY <- cbind(YY,
              ifelse(j == 1, "Anchored yields only DNS  ", "Anchored macro augmented DNS"))
  colnames(YY) <- c("Dates", "h", "tau", "rm", "mm")
  YY$Dates <- dat[-c(1:47)]
  roll_mod_info[[j]] <- bind_rows(XX, YY)
}

roll_mod_info_2 <- list()
for (j in 1:2) {
  XX <- roll_rel_msfe[[1]][[j]][, hh, ggmat]
  XX <- melt(XX)
  XX <- cbind(XX, ifelse(j == 1, "Yields only DNS", "Macro augmented DNS"))
  colnames(XX) <- c("Dates", "h", "tau", "rm", "mm")
  XX$Dates <- dat[-c(1:47)]
  YY <- roll_rel_msfe[[3]][[j]][, hh, ggmat]
  YY <- melt(YY)
  YY <- cbind(YY,
              ifelse(j == 1, "Anchored yields only DNS (3m adn 6m)", "Anchored macro augmented DNS (3m and 6m)"))
  colnames(YY) <- c("Dates", "h", "tau", "rm", "mm")
  YY$Dates <- dat[-c(1:47)]
  roll_mod_info_2[[j]] <- bind_rows(XX, YY)
}


title_graphs = list(a = "(a) Yields only DNS",
                    b = "(b) Macro augmented DNS")

```


```{r calculate_mae, echo=FALSE}
## dns_model is a list each is a model
maturities <- c(seq(3,24,3), 30, seq(36, 120, 12))
cnames <- paste0('$\\tau=', maturities, "$")
rnames <- paste0('$h=', 1:12, "$")

mae_ss <- map(dns_model, ~MAE(.x$error, colnames = cnames, rownames = rnames))
mae_rw <- map(dns_model, ~MAE(.x$error_rw, colnames = cnames, rownames = rnames))
mae_tl <- map(dns_model, ~MAE(.x$error_tilt, colnames = cnames, rownames = rnames))
mae_tl_2 <- map(dns_model, ~MAE(.x$error_tilt_2, colnames = cnames, rownames = rnames))

mae_cf <- map(dns_model, ~MAE(.x$error_cond, colnames = cnames, rownames = rnames))

rel_mae <- list(
  ss_ss = map2(mae_ss, mae_rw, ~ .x / .y),
  ss_tl = map2(mae_tl, mae_rw, ~ .x / .y),
  ss_tl_2 = map2(mae_tl_2, mae_rw, ~ .x / .y),
  ss_cf = map2(mae_cf, mae_rw, ~ .x / .y)
  )

mae_ss <- map(dns_model, ~MAE(.x$error[1:108, , ],
                                colnames = cnames, rownames = rnames))
mae_rw <- map(dns_model, ~MAE(.x$error_rw[1:108, , ],
                                colnames = cnames, rownames = rnames))
mae_tl <- map(dns_model, ~MAE(.x$error_tilt[1:108, , ],
                                colnames = cnames, rownames = rnames))

mae_tl_2 <- map(dns_model, ~MAE(.x$error_tilt_2[1:108, , ],
                                colnames = cnames, rownames = rnames))

mae_cf <- map(dns_model, ~MAE(.x$error_cond[1:108, , ],
                                colnames = cnames, rownames = rnames))

rel_mae_up2008 <- list(
  ss_ss = map2(mae_ss, mae_rw, ~ .x / .y),
  ss_tl = map2(mae_tl, mae_rw, ~ .x / .y),
  ss_tl_2 = map2(mae_tl_2, mae_rw, ~ .x / .y),
  ss_cf = map2(mae_cf, mae_rw, ~ .x / .y)
  )


mae_ss <- map(dns_model, ~MAE(.x$error[49:142, , ],
                                colnames = cnames, rownames = rnames))
mae_rw <- map(dns_model, ~MAE(.x$error_rw[49:142, , ],
                                colnames = cnames, rownames = rnames))
mae_tl <- map(dns_model, ~MAE(.x$error_tilt[49:142, , ],
                                colnames = cnames, rownames = rnames))

mae_tl_2 <- map(dns_model, ~MAE(.x$error_tilt_2[49:142, , ],
                                colnames = cnames, rownames = rnames))

mae_cf <- map(dns_model, ~MAE(.x$error_cond[49:142, , ],
                                colnames = cnames, rownames = rnames))

rel_mae_from2005 <- list(
  ss_ss = map2(mae_ss, mae_rw, ~ .x / .y),
  ss_tl = map2(mae_tl, mae_rw, ~ .x / .y),
  ss_tl_2 = map2(mae_tl_2, mae_rw, ~ .x / .y),
  ss_cf = map2(mae_cf, mae_rw, ~ .x / .y)
  )


mae_ss <- map(dns_model, ~RollMAE(.x$error, n = 48,
                                    colnames = cnames, rownames = rnames))
mae_rw <- map(dns_model, ~RollMAE(.x$error_rw, n = 48,
                                    colnames = cnames, rownames = rnames))
mae_tl <- map(dns_model, ~RollMAE(.x$error_tilt, n = 48,
                                    colnames = cnames, rownames = rnames))
mae_tl_2 <- map(dns_model, ~RollMAE(.x$error_tilt_2, n = 48,
                                    colnames = cnames, rownames = rnames))
mae_cf <- map(dns_model, ~RollMAE(.x$error_cond, n = 48,
                                    colnames = cnames, rownames = rnames))

roll_rel_mae <- list(
  ss_ss = map2(mae_ss, mae_rw, ~ .x / .y),
  ss_tl = map2(mae_tl, mae_rw, ~ .x / .y),
  ss_tl_2 = map2(mae_tl_2, mae_rw, ~ .x / .y),
  ss_cf = map2(mae_cf, mae_rw, ~ .x / .y)
  )


## Create rollin matrix to be plotted using ggplot
ggmat <- c(3, 6, 9, 12, 17)
hh <- c(3, 6, 9, 12)
dat <- seq(ymd("2000-01-2"), ymd("2012-01-2"), by = "months")
roll_mae_info <- list()
for (j in 1:2) {
  XX <- roll_rel_mae[[1]][[j]][, hh, ggmat]
  XX <- melt(XX)
  XX <- cbind(XX, ifelse(j == 1, "Yields only DNS", "Macro augmented DNS"))
  colnames(XX) <- c("Dates", "h", "tau", "rm", "mm")
  XX$Dates <- dat[-c(1:47)]
  YY <- roll_rel_mae[[2]][[j]][, hh, ggmat]
  YY <- melt(YY)
  YY <- cbind(YY,
              ifelse(j == 1, "Anchored yields only DNS  ", "Anchored macro augmented DNS"))
  colnames(YY) <- c("Dates", "h", "tau", "rm", "mm")
  YY$Dates <- dat[-c(1:47)]
  roll_mae_info[[j]] <- bind_rows(XX, YY)
}

roll_mae_info_2 <- list()
for (j in 1:2) {
  XX <- roll_rel_msfe[[1]][[j]][, hh, ggmat]
  XX <- melt(XX)
  XX <- cbind(XX, ifelse(j == 1, "Yields only DNS", "Macro augmented DNS"))
  colnames(XX) <- c("Dates", "h", "tau", "rm", "mm")
  XX$Dates <- dat[-c(1:47)]
  YY <- roll_rel_mae[[3]][[j]][, hh, ggmat]
  YY <- melt(YY)
  YY <- cbind(YY,
              ifelse(j == 1, "Anchored yields only DNS (3m adn 6m)", "Anchored macro augmented DNS (3m and 6m)"))
  colnames(YY) <- c("Dates", "h", "tau", "rm", "mm")
  YY$Dates <- dat[-c(1:47)]
  roll_mae_info_2[[j]] <- bind_rows(XX, YY)
}


title_graphs = list(a = "(a) Yields only DNS",
                    b = "(b) Macro augmented DNS")

```


```{r relative_mse_dns, echo=FALSE, fig.cap="Relative MSFE of the forecasts based on the  yields only and macro augmented DNS models against the random walk forecasts"}

par(
  mfrow = c(1, 2),
  mar = c(3, 3, 2, 1),
  mgp = c(2, .7, 0),
  tck = -.01,
  cex.lab = 0.75,
  cex.axis = 0.75
  )

walk2(rel_msfe[[1]], title_graphs, ~ plot_msfe(.x, .y, color = FALSE))
```

```{r relative_msfe_dns_tilted, fig.cap="Relative MSFE of the tilted forecasts based on the  yields only and macro augmented DNS models against the random walk forecasts"}

par(
  mfrow = c(1, 2),
  mar = c(3, 3, 2, 1),
  mgp = c(2, .7, 0),
  tck = -.01,
  cex.lab = 0.75,
  cex.axis = 0.75
  )

walk2(rel_msfe[[2]], title_graphs, ~plot_msfe(.x, .y, color = FALSE))

```


```{r rolling_base, fig.width = 6.2, fig.height = 7, fig.cap="Rolling relative MSFE of the forecasts based on the  yields only and macro augmented DNS models against the random walk forecasts"}
for (j in 1) {
  print(
    ggplot(roll_mod_info[[j]], aes(y = rm, x = Dates)) +
      geom_line(aes(col = mm, lty = mm), size = .6) +
      facet_grid(tau ~ h, switch = "y") +
      theme_bw(base_size = 8) + ylab("Mean squared forecasting error") + xlab("Time") +
      theme(
        panel.border = element_rect(linetype = "solid", colour = "black", size = 1.2),
        legend.position = "bottom",
        legend.title = element_blank(),
        strip.background = element_rect(fill = "#FFFFFF", color = "#FFFFFF", size = 1),
        strip.text = element_text(size = 7)) +
      scale_color_manual(values = c("#000000", "#000000")) +
      scale_linetype_manual(values = c("solid", "dashed"))

      )
}
```

```{r rolling_macro, fig.width = 6.2, fig.height = 7, fig.cap="Rolling relative MSFE of the tilted forecasts based on the  yields only and macro augmented DNS models against the random walk forecasts"}
for (j in 2) {
  print(
    ggplot(roll_mod_info[[j]], aes(y = rm, x = Dates)) +
      geom_line(aes(col = mm, lty = mm), size = .6) +
      facet_grid(tau ~ h, switch = "y") +
      theme_bw(base_size = 8) + ylab("Mean squared forecasting error") + xlab("Time") +
      theme(
        panel.border = element_rect(linetype = "solid", colour = "black", size = 1.2),        legend.position = "bottom",
        legend.title = element_blank(),
        strip.background = element_rect(fill = "#FFFFFF", color = "#FFFFFF", size = 1),
        strip.text = element_text(size = 6)) +
      scale_color_manual(values = c("#000000", "#000000")) +
      scale_linetype_manual(values = c("solid", "dashed"))
      )
}
```


```{r rolling_macro_2, fig.width = 6.2, fig.height = 7, fig.cap="Rolling relative MSFE of the tilted forecasts based on the  yields only and macro augmented DNS models against the random walk forecasts. Both the 3-month and the 6-month BC forecasts."}
for (j in 2) {
  print(
    ggplot(roll_mod_info_2[[j]], aes(y = rm, x = Dates)) +
      geom_line(aes(col = mm, lty = mm), size = .6) +
      facet_grid(tau ~ h, switch = "y") +
      theme_bw(base_size = 8) + ylab("Mean squared forecasting error") + xlab("Time") +
      theme(
        panel.border = element_rect(linetype = "solid", colour = "black", size = 1.2),        legend.position = "bottom",
        legend.title = element_blank(),
        strip.background = element_rect(fill = "#FFFFFF", color = "#FFFFFF", size = 1),
        strip.text = element_text(size = 6)) +
      scale_color_manual(values = c("#000000", "#000000")) +
      scale_linetype_manual(values = c("solid", "dashed"))
      )
}
```


```{r rolling_macro_mae, fig.width = 6.2, fig.height = 7, fig.cap="Rolling relative MAFE of the tilted forecasts based on the  yields only and macro augmented DNS models against the random walk forecasts"}
for (j in 2) {
  print(
    ggplot(roll_mae_info[[j]], aes(y = rm, x = Dates)) +
      geom_line(aes(col = mm, lty = mm), size = .6) +
      facet_grid(tau ~ h, switch = "y") +
      theme_bw(base_size = 8) + ylab("Mean absolute forecasting error") + xlab("Time") +
      theme(
        panel.border = element_rect(linetype = "solid", colour = "black", size = 1.2),        legend.position = "bottom",
        legend.title = element_blank(),
        strip.background = element_rect(fill = "#FFFFFF", color = "#FFFFFF", size = 1),
        strip.text = element_text(size = 6)) +
      scale_color_manual(values = c("#000000", "#000000")) +
      scale_linetype_manual(values = c("solid", "dashed"))
      )
}
```


```{r rolling_macro_mae_2, fig.width = 6.2, fig.height = 7, fig.cap="Rolling relative MAFE of the tilted forecasts based on the  yields only and macro augmented DNS models against the random walk forecasts. Both the 3-month and the 6-month BC forecasts."}
for (j in 2) {
  print(
    ggplot(roll_mod_info_2[[j]], aes(y = rm, x = Dates)) +
      geom_line(aes(col = mm, lty = mm), size = .6) +
      facet_grid(tau ~ h, switch = "y") +
      theme_bw(base_size = 8) + ylab("Mean absolute forecasting error") + xlab("Time") +
      theme(
        panel.border = element_rect(linetype = "solid", colour = "black", size = 1.2),        legend.position = "bottom",
        legend.title = element_blank(),
        strip.background = element_rect(fill = "#FFFFFF", color = "#FFFFFF", size = 1),
        strip.text = element_text(size = 6)) +
      scale_color_manual(values = c("#000000", "#000000")) +
      scale_linetype_manual(values = c("solid", "dashed"))
      )
}
```



```{r encompassing, echo=FALSE, fig.height = 4.2, fig.cap="Notes: The figure reports the sequence of test statistics for the time-varying encompassing test described in Section 3.2, testing the null hypothesis that the BC forecast encompasses the DNS forecast, against the alternative hypothesis that it does not. The null hypothesis is rejected when the sequence of test statistics crosses the horizontal solid line, which represents the critical value (which equals 2.62 for test statistics computed over an estimation window that uses $40\\%$ of the out-of-sample observations and for a $5\\%$ significance level)."}


off <- which(index(yield_data) == "Jan 2000")
bc_error <- list()
mat_ <- c(1, 2, 4, 8, 12, 17)
excl <- -c(1:(off - 1))

for (mm in c(3, 6, 12, 24, 60, 120)) {
  mat <- paste0("m", mm)
  mu  <- bc_set[[mat]][[2]]
  idx <- match(mm, maturities)
  for (hh in c(3, 6, 9, 12)) {
    bc_error[[mat]][[paste0("h=", hh)]] <-
      yield_data[excl, idx] - mu[[hh / 3]][excl]
  }
}
msfe_bc <- matrix(0, 6, 4)
for (mm in 1:6) {
  msfe_bc[mm, ] <- sapply(bc_error[[mm]], MSFE)
}

rel_msfe_bc_rw <-
  msfe_bc / t(MSFE(dns_model[[1]]$error_rw)[c(3, 6, 9, 12), c(1, 2, 4, 8, 12, 17)])

colnames(rel_msfe_bc_rw) <-
  paste0("$h=", c(3, 6, 9, 12), "$")
rownames(rel_msfe_bc_rw) <-
  paste0("$\\tau=", c(3, 6, 12, 24, 60, 120), "$")

plot_encompassing_2()


#x <- dns_model[[1]]$error; Sigma <-  dns_model[[1]]$Sigma
```
\clearpage


```{r encompassing_36, echo=FALSE, fig.height = 4.2, fig.cap="Notes: The figure reports the sequence of test statistics for the time-varying encompassing test described in Section 3.2, testing the null hypothesis that the BC forecast encompasses the DNS forecast, against the alternative hypothesis that it does not. The null hypothesis is rejected when the sequence of test statistics crosses the horizontal solid line, which represents the critical value (which equals 2.62 for test statistics computed over an estimation window that uses $40\\%$ of the out-of-sample observations and for a $5\\%$ significance level)."}


plot_encompassing_3()


#x <- dns_model[[1]]$error; Sigma <-  dns_model[[1]]$Sigma
```
\clearpage




```{r relative_mse_base_full,  echo = FALSE, results='asis'}

## (One sided)
## The null hypothesis is
## H0: MSE(e1) ≤ MSE(e2)
## H1: MSE(e1) ≥ MSE(e2)  

putstars <- function(x) {
  stars <- matrix("", nrow(x), ncol(x))
  stars[x <= 0.01] <- "***"
  stars[0.01 < x & x <= 0.05] <- "**"
  stars[0.05 < x & x <= 0.1] <- "*"
  stars
}

tbl_list <- list()

for (jj in 1:2) {
  tbl_1 <- matrix(0, 17,4)
  for (j in seq(3, 12, 3)) {
    tbl_1[,j/3] <- apply(dns_model[[jj]]$error_tilt[,j,], 2, MSFE)/
      apply(dns_model[[jj]]$error[,j,], 2, MSFE)
  }

  tbl_2 <- matrix(0, 17,4)
  for (j in seq(3, 12, 3)) {
    tbl_2[,j/3] <- apply(dns_model[[jj]]$error_tilt[,j,], 2, MSFE)/
      apply(dns_model[[jj]]$error_rw[,j,], 2, MSFE)
  }

  pval_1 <- pval_2 <- NULL
  for (j in seq(3, 12, 3)) {
    e1 <- as.data.frame(dns_model[[jj]]$error_rw[,j,])
    e2 <- as.data.frame(dns_model[[jj]]$error_tilt[,j,])
    e3 <- as.data.frame(dns_model[[jj]]$error[,j,])
    mdl_1 <- map2(e1, e2, dm_test, alternative = "greater", h = j)
    mdl_2 <- map2(e2, e3, dm_test, alternative = "greater", h = j)
    pval_1 <- cbind(pval_1, map_dbl(mdl_1, ~.x$p.value))
    pval_2 <- cbind(pval_2, map_dbl(mdl_1, ~.x$p.value))
  }

  stars_1 <- putstars(pval_1)
  stars_2 <- putstars(pval_2)
  tbl_1 <- format(tbl_1, digits = 2)
  for (j in 1:nrow(tbl_1)) {
    for (k in 1:ncol(tbl_1))
      tbl_1[j,k] <- paste0("$", tbl_1[j,k], "^{", stars_1[j,k], "}$")
  }

  tbl_2 <- format(tbl_2, digits = 2)
  for (j in 1:nrow(tbl_2)) {
    for (k in 1:ncol(tbl_2))
      tbl_2[j,k] <- paste0("$", tbl_2[j,k], "^{", stars_2[j,k], "}$")
  }

  out <- cbind(tbl_1, tbl_2)
  rownames(out) <- paste(maturities)
  colnames(out) <- rep(c('h=3', 'h=6', 'h=9', 'h=12'),2)

  tbl_list[[jj]] <- list(out,
                         description = dns_model[[jj]]$description)
}

options(latexcmd = 'pdflatex', xdvicmd = 'open',  dviExtension = 'pdf')

for (jj in c(1, 2)) {
  latex.def(
    tbl_list[[jj]][[1]],
    dcolumn = TRUE,
    file = "",
    title = "",
    rowlabel = "Maturity",
    cgroup = c("Anchored vs DNS", "Anchored vs RW"),
    n.cgroup = c(4, 4),
    dec = 2,
    align = "crrrrcrrrr",
    numeric.dollar = FALSE,
    size = "footnotesize",
    caption = paste0(
      paste0(
        "\\textbf{Relative MSFEs of anchored forecasts. Baseline model:",
        tbl_list[[jj]]$description
      ),
      " Sample: 2000:01 - 2011:12.} The table reports the ratios of MSFE for the models considered. The asterisk indicates significance  according to the Diebold and Mariano (1995) test of equal accuracy against the alternative that the anchored forecast is more accurate ('***' at the 1\\%, '**' at the 5\\%, and '*' at the 1\\%). The Diebold and Mariano test was implemented using an HAC estimator with $h-1$ truncation parameter."
    ),
    caption.loc = "bottom",
    nomargin = FALSE,
    comment = FALSE
  )

  nms <-
    c("tables/relative_mse_base.tex",
      "tables/relative_mse_macro.tex")

  w <- latex.def(
      tbl_list[[jj]][[1]],
      dcolumn = TRUE,
      file = nms[jj],
      title = "",
      table.env = FALSE,
      rowlabel = "Maturity",
      cgroup = c("Anchored vs DNS", "Anchored vs RW"),
      n.cgroup = c(4, 4),
      dec = 2,
      align = "crrrrcrrrr",
      numeric.dollar = FALSE,
      size = "footnotesize")
}

```


```{r relative_mse_base_full_2,  echo = FALSE, results='asis'}

## (One sided)
## The null hypothesis is
## H0: MSE(e1) ≤ MSE(e2)
## H1: MSE(e1) ≥ MSE(e2)  

putstars <- function(x) {
  stars <- matrix("", nrow(x), ncol(x))
  stars[x <= 0.01] <- "***"
  stars[0.01 < x & x <= 0.05] <- "**"
  stars[0.05 < x & x <= 0.1] <- "*"
  stars
}

tbl_list <- list()

for (jj in 1:2) {
  tbl_1 <- matrix(0, 17,4)
  for (j in seq(3, 12, 3)) {
    tbl_1[,j/3] <- apply(dns_model[[jj]]$error_tilt_2[,j,], 2, MSFE)/
      apply(dns_model[[jj]]$error[,j,], 2, MSFE)
  }

  tbl_2 <- matrix(0, 17,4)
  for (j in seq(3, 12, 3)) {
    tbl_2[,j/3] <- apply(dns_model[[jj]]$error_tilt_2[,j,], 2, MSFE)/
      apply(dns_model[[jj]]$error_rw[,j,], 2, MSFE)
  }

  pval_1 <- pval_2 <- NULL
  for (j in seq(3, 12, 3)) {
    e1 <- as.data.frame(dns_model[[jj]]$error_rw[,j,])
    e2 <- as.data.frame(dns_model[[jj]]$error_tilt_2[,j,])
    e3 <- as.data.frame(dns_model[[jj]]$error[,j,])
    mdl_1 <- map2(e1, e2, dm_test, alternative = "greater", h = j)
    mdl_2 <- map2(e2, e3, dm_test, alternative = "greater", h = j)
    pval_1 <- cbind(pval_1, map_dbl(mdl_1, ~.x$p.value))
    pval_2 <- cbind(pval_2, map_dbl(mdl_1, ~.x$p.value))
  }

  stars_1 <- putstars(pval_1)
  stars_2 <- putstars(pval_2)
  tbl_1 <- format(tbl_1, digits = 2)
  for (j in 1:nrow(tbl_1)) {
    for (k in 1:ncol(tbl_1))
      tbl_1[j,k] <- paste0("$", tbl_1[j,k], "^{", stars_1[j,k], "}$")
  }

  tbl_2 <- format(tbl_2, digits = 2)
  for (j in 1:nrow(tbl_2)) {
    for (k in 1:ncol(tbl_2))
      tbl_2[j,k] <- paste0("$", tbl_2[j,k], "^{", stars_2[j,k], "}$")
  }

  out <- cbind(tbl_1, tbl_2)
  rownames(out) <- paste(maturities)
  colnames(out) <- rep(c('h=3', 'h=6', 'h=9', 'h=12'),2)

  tbl_list[[jj]] <- list(out,
                         description = dns_model[[jj]]$description)
}

options(latexcmd = 'pdflatex', xdvicmd = 'open',  dviExtension = 'pdf')

for (jj in c(1, 2)) {
  latex.def(
    tbl_list[[jj]][[1]],
    dcolumn = TRUE,
    file = "",
    title = "",
    rowlabel = "Maturity",
    cgroup = c("Anchored vs DNS", "Anchored vs RW"),
    n.cgroup = c(4, 4),
    dec = 2,
    align = "crrrrcrrrr",
    numeric.dollar = FALSE,
    size = "footnotesize",
    caption = paste0(
      paste0(
        "\\textbf{Relative MSFEs of anchored forecasts using both 3m and 6m BC forecasts. Baseline model:",
        tbl_list[[jj]]$description
      ),
      " Sample: 2000:01 - 2011:12.} The table reports the ratios of MSFE for the models considered. The asterisk indicates significance  according to the Diebold and Mariano (1995) test of equal accuracy against the alternative that the anchored forecast is more accurate ('***' at the 1\\%, '**' at the 5\\%, and '*' at the 1\\%). The Diebold and Mariano test was implemented using an HAC estimator with $h-1$ truncation parameter."
    ),
    caption.loc = "bottom",
    nomargin = FALSE,
    comment = FALSE
  )

  nms <-
    c("tables/relative_mse_base_2.tex",
      "tables/relative_mse_macro_2.tex")

  w <- latex.def(
      tbl_list[[jj]][[1]],
      dcolumn = TRUE,
      file = nms[jj],
      title = "",
      table.env = FALSE,
      rowlabel = "Maturity",
      cgroup = c("Anchored vs DNS", "Anchored vs RW"),
      n.cgroup = c(4, 4),
      dec = 2,
      align = "crrrrcrrrr",
      numeric.dollar = FALSE,
      size = "footnotesize")
}

```

\clearpage
```{r relative_mse_base_upto2008, echo = FALSE, results='asis'}
## (One sided)
## The null hypothesis is
## H0: MSE(e1) ≤ MSE(e2)
## H1: MSE(e1) ≥ MSE(e2)  

putstars <- function(x) {
  stars <- matrix("", nrow(x), ncol(x))
  stars[x <= 0.01] <- "***"
  stars[0.01 < x & x <= 0.05] <- "**"
  stars[0.05 < x & x <= 0.1] <- "*"
  stars
  }
idxx <- 1:108
tbl_list <- list()

for (jj in 1:2) {
  tbl_1 <- matrix(0, 17, 4)
  for (j in seq(3, 12, 3)) {
    tbl_1[,j/3] <- apply(dns_model[[jj]]$error_tilt[idxx, j, ], 2, MSFE)/
      apply(dns_model[[jj]]$error[idxx, j,], 2, MSFE)
  }

  tbl_2 <- matrix(0, 17,4)
  for (j in seq(3, 12, 3)) {
    tbl_2[,j/3] <- apply(dns_model[[jj]]$error_tilt[idxx, j,], 2, MSFE)/
      apply(dns_model[[jj]]$error_rw[idxx, j,], 2, MSFE)
  }

  pval_1 <- pval_2 <- NULL
  for (j in seq(3, 12, 3)) {
    e1 <- as.data.frame(dns_model[[jj]]$error_rw[idxx, j,])
    e2 <- as.data.frame(dns_model[[jj]]$error_tilt[idxx, j,])
    e3 <- as.data.frame(dns_model[[jj]]$error[idxx,j,])
    mdl_1 <- map2(e1, e2, dm_test, alternative = "greater", h = j)
    mdl_2 <- map2(e2, e3, dm_test, alternative = "greater", h = j)
    pval_1 <- cbind(pval_1, map_dbl(mdl_1, ~.x$p.value))
    pval_2 <- cbind(pval_2, map_dbl(mdl_1, ~.x$p.value))
  }

  stars_1 <- putstars(pval_1)
  stars_2 <- putstars(pval_2)
  tbl_1 <- format(tbl_1, digits = 2)
  for (j in 1:nrow(tbl_1)) {
    for (k in 1:ncol(tbl_1))
      tbl_1[j,k] <- paste0("$", tbl_1[j,k], "^{", stars_1[j,k], "}$")
  }

  tbl_2 <- format(tbl_2, digits = 2)
  for (j in 1:nrow(tbl_2)) {
    for (k in 1:ncol(tbl_2))
      tbl_2[j,k] <- paste0("$", tbl_2[j,k], "^{", stars_2[j,k], "}$")
  }

  out <- cbind(tbl_1, tbl_2)
  rownames(out) <- paste(maturities)
  colnames(out) <- rep(c('h=3', 'h=6', 'h=9', 'h=12'),2)

  tbl_list[[jj]] <- list(out,
                         description = dns_model[[jj]]$description)
}



for (jj in 1:2) {

  latex.def(
    tbl_list[[jj]][[1]],
    dcolumn = TRUE,
    file = "",
    title = "",
    rowlabel = "Maturity",
    cgroup = c("Anchored vs DNS", "Anchored vs RW"),
    n.cgroup = c(4, 4),
    dec = 2,
    align = "crrrrcrrrr",
    numeric.dollar = FALSE,
    size = "footnotesize",
    caption = paste0(
    paste0(
    "\\textbf{Relative MSFEs of anchored forecasts. DNS model:",
    tbl_list[[jj]]$description
    ),
    " Sample: 2000:01 - 2008:12.} The table reports the ratios of MSFE for the models considered. The asterisk indicates significance  according to the Diebold and Mariano (1995) test of equal accuracy against the alternative that the anchored forecast is more accurate ('***' at the 1\\%, '**' at the 5\\%, and '*' at the 1\\%). The Diebold and Mariano test was implemented using an HAC estimator with $h-1$ truncation parameter."
    ),
    caption.loc = "bottom",
    nomargin = FALSE,
    comment = FALSE
    )

   nms <-
    c("tables/relative_mse_base_upto2008.tex",
      "tables/relative_mse_macro_upto2008.tex")

  w <- latex.def(
    tbl_list[[jj]][[1]],
    dcolumn = TRUE,
    file = nms[jj],
    title = "",
    table.env = FALSE,
    rowlabel = "Maturity",
    cgroup = c("Anchored vs DNS", "Anchored vs RW"),
    n.cgroup = c(4, 4),
    dec = 2,
    align = "crrrrcrrrr",
    numeric.dollar = FALSE,
    size = "footnotesize"
    )

}  
```

```{r relative_mse_base_upto2008_2, echo = FALSE, results='asis'}
## (One sided)
## The null hypothesis is
## H0: MSE(e1) ≤ MSE(e2)
## H1: MSE(e1) ≥ MSE(e2)  

putstars <- function(x) {
  stars <- matrix("", nrow(x), ncol(x))
  stars[x <= 0.01] <- "***"
  stars[0.01 < x & x <= 0.05] <- "**"
  stars[0.05 < x & x <= 0.1] <- "*"
  stars
  }
idxx <- 1:108
tbl_list <- list()

for (jj in 1:2) {
  tbl_1 <- matrix(0, 17, 4)
  for (j in seq(3, 12, 3)) {
    tbl_1[,j/3] <- apply(dns_model[[jj]]$error_tilt_2[idxx, j, ], 2, MSFE)/
      apply(dns_model[[jj]]$error[idxx, j,], 2, MSFE)
  }

  tbl_2 <- matrix(0, 17,4)
  for (j in seq(3, 12, 3)) {
    tbl_2[,j/3] <- apply(dns_model[[jj]]$error_tilt_2[idxx, j,], 2, MSFE)/
      apply(dns_model[[jj]]$error_rw[idxx, j,], 2, MSFE)
  }

  pval_1 <- pval_2 <- NULL
  for (j in seq(3, 12, 3)) {
    e1 <- as.data.frame(dns_model[[jj]]$error_rw[idxx, j,])
    e2 <- as.data.frame(dns_model[[jj]]$error_tilt_2[idxx, j,])
    e3 <- as.data.frame(dns_model[[jj]]$error[idxx,j,])
    mdl_1 <- map2(e1, e2, dm_test, alternative = "greater", h = j)
    mdl_2 <- map2(e2, e3, dm_test, alternative = "greater", h = j)
    pval_1 <- cbind(pval_1, map_dbl(mdl_1, ~.x$p.value))
    pval_2 <- cbind(pval_2, map_dbl(mdl_1, ~.x$p.value))
  }

  stars_1 <- putstars(pval_1)
  stars_2 <- putstars(pval_2)
  tbl_1 <- format(tbl_1, digits = 2)
  for (j in 1:nrow(tbl_1)) {
    for (k in 1:ncol(tbl_1))
      tbl_1[j,k] <- paste0("$", tbl_1[j,k], "^{", stars_1[j,k], "}$")
  }

  tbl_2 <- format(tbl_2, digits = 2)
  for (j in 1:nrow(tbl_2)) {
    for (k in 1:ncol(tbl_2))
      tbl_2[j,k] <- paste0("$", tbl_2[j,k], "^{", stars_2[j,k], "}$")
  }

  out <- cbind(tbl_1, tbl_2)
  rownames(out) <- paste(maturities)
  colnames(out) <- rep(c('h=3', 'h=6', 'h=9', 'h=12'),2)

  tbl_list[[jj]] <- list(out,
                         description = dns_model[[jj]]$description)
}



for (jj in 1:2) {

  latex.def(
    tbl_list[[jj]][[1]],
    dcolumn = TRUE,
    file = "",
    title = "",
    rowlabel = "Maturity",
    cgroup = c("Anchored vs DNS", "Anchored vs RW"),
    n.cgroup = c(4, 4),
    dec = 2,
    align = "crrrrcrrrr",
    numeric.dollar = FALSE,
    size = "footnotesize",
    caption = paste0(
    paste0(
    "\\textbf{Relative MSFEs of anchored forecasts using 3m and 6m BC survey forecasts. DNS model:",
    tbl_list[[jj]]$description
    ),
    " Sample: 2000:01 - 2008:12.} The table reports the ratios of MSFE for the models considered. The asterisk indicates significance  according to the Diebold and Mariano (1995) test of equal accuracy against the alternative that the anchored forecast is more accurate ('***' at the 1\\%, '**' at the 5\\%, and '*' at the 1\\%). The Diebold and Mariano test was implemented using an HAC estimator with $h-1$ truncation parameter."
    ),
    caption.loc = "bottom",
    nomargin = FALSE,
    comment = FALSE
    )

   nms <-
    c("tables/relative_mse_base_upto2008_2.tex",
      "tables/relative_mse_macro_upto2008_2.tex")

  w <- latex.def(
    tbl_list[[jj]][[1]],
    dcolumn = TRUE,
    file = nms[jj],
    title = "",
    table.env = FALSE,
    rowlabel = "Maturity",
    cgroup = c("Anchored vs DNS", "Anchored vs RW"),
    n.cgroup = c(4, 4),
    dec = 2,
    align = "crrrrcrrrr",
    numeric.dollar = FALSE,
    size = "footnotesize"
    )

}  
```

```{r relative_mse_base_full_cond, echo = FALSE, results='asis'}

## (One sided)
## The null hypothesis is
## H0: MSE(e1) ≤ MSE(e2)
## H1: MSE(e1) ≥ MSE(e2)  

putstars <- function(x) {
  stars <- matrix("", nrow(x), ncol(x))
  stars[x <= 0.01] <- "***"
  stars[0.01 < x & x <= 0.05] <- "**"
  stars[0.05 < x & x <= 0.1] <- "*"
  stars
}

tbl_list <- list()

for (jj in 1:2) {
  tbl_1 <- matrix(0, 17,4)
  for (j in seq(3, 12, 3)) {
    tbl_1[,j/3] <- apply(dns_model[[jj]]$error_cond[,j,], 2, MSFE)/
      apply(dns_model[[jj]]$error[,j,], 2, MSFE)
  }

  tbl_2 <- matrix(0, 17,4)
  for (j in seq(3, 12, 3)) {
    tbl_2[,j/3] <- apply(dns_model[[jj]]$error_cond[,j,], 2, MSFE)/
      apply(dns_model[[jj]]$error_rw[,j,], 2, MSFE)
  }

  pval_1 <- pval_2 <- NULL
  for (j in seq(3, 12, 3)) {
    e1 <- as.data.frame(dns_model[[jj]]$error_rw[,j,])
    e2 <- as.data.frame(dns_model[[jj]]$error_cond[,j,])
    e3 <- as.data.frame(dns_model[[jj]]$error[,j,])
    mdl_1 <- map2(e1, e2, dm_test, alternative = "greater", h = j)
    mdl_2 <- map2(e2, e3, dm_test, alternative = "greater", h = j)
    pval_1 <- cbind(pval_1, map_dbl(mdl_1, ~.x$p.value))
    pval_2 <- cbind(pval_2, map_dbl(mdl_1, ~.x$p.value))
  }

  stars_1 <- putstars(pval_1)
  stars_2 <- putstars(pval_2)
  tbl_1 <- format(tbl_1, digits = 2)
  for (j in 1:nrow(tbl_1)) {
    for (k in 1:ncol(tbl_1))
      tbl_1[j,k] <- paste0("$", tbl_1[j,k], "^{", stars_1[j,k], "}$")
  }

  tbl_2 <- format(tbl_2, digits = 2)
  for (j in 1:nrow(tbl_2)) {
    for (k in 1:ncol(tbl_2))
      tbl_2[j,k] <- paste0("$", tbl_2[j,k], "^{", stars_2[j,k], "}$")
  }

  out <- cbind(tbl_1, tbl_2)
  rownames(out) <- paste(maturities)
  colnames(out) <- rep(c('h=3', 'h=6', 'h=9', 'h=12'),2)

  tbl_list[[jj]] <- list(out,
                         description = dns_model[[jj]]$description)
}



for (jj in 1:2) {

  latex.def(
    tbl_list[[jj]][[1]],
    dcolumn = TRUE,
    file = "",
    title = "",
    rowlabel = "Maturity",
    cgroup = c("Conditional forecast vs DNS", "Conditional vs RW"),
    n.cgroup = c(4, 4),
    dec = 2,
    align = "crrrrcrrrr",
    numeric.dollar = FALSE,
    size = "footnotesize",
    caption = paste0(
    paste0(
    "\\textbf{Relative MSFEs of conditional forecasts. Baseline model:",
    tbl_list[[jj]]$description
    ),
    " Sample: 2000:01 - 2011:12.} The table reports the ratios of MSFE for the models considered. The asterisk indicates significance  according to the Diebold and Mariano (1995) test of equal accuracy against the alternative that the anchored forecast is more accurate ('***' at the 1\\%, '**' at the 5\\%, and '*' at the 1\\%). The Diebold and Mariano test was implemented using an HAC estimator with $h-1$ truncation parameter."
    ),
    caption.loc = "bottom",
    nomargin = FALSE,
    comment = FALSE
    )


  nms <- c("tables/relative_mse_conditional_base.tex",
           "tables/relative_mse_conditional_macro.tex")

  w <- latex.def(
    tbl_list[[jj]][[1]],
    dcolumn = TRUE,
    file = nms[jj],
    title = "",
    table.env = FALSE,
    rowlabel = "Maturity",
    cgroup = c("Conditional vs DNS", "Conditional vs RW"),
    n.cgroup = c(4, 4),
    dec = 2,
    align = "crrrrcrrrr",
    numeric.dollar = FALSE,
    size = "footnotesize"
    )
  }

```



```{r fomc, fig.cap="The figure reports the 1- to 4-quarter-ahead forecasts of the 3-month yield given by the DNS model and the BC survey before and after the FOMC Statement of August 9, 2011."}

## 145 - December
## 144 - November
## 143 - October
## 142 - September
## 141 - August
## 140 - July
## 139 - June
## 138 - May
## 137 - April
## 136 - March
## 135 - Februry
## 134 - January

par(
  mfcol = c(1,1),
  mar = c(3, 3, 2, 1),
  mgp = c(2, .7, 0),
  tck = -.01,
  cex.lab = 0.75,
  cex.axis = 0.75
)

didx <- index(yield_data)
idxx <- which(didx == "Jun 2011") - 180
xc <- 1:29
yy <- c(dns_model[[1]]$fcast[idxx + 1, c(3, 6, 9, 12),1],
        NA,
        dns_model[[1]]$fcast[idxx + 2, c(3, 6, 9, 12), 1],
        NA,
        dns_model[[1]]$fcast[idxx + 3, c(3, 6, 9, 12), 1],
        NA,
        dns_model[[1]]$fcast[idxx + 4, c(3, 6, 9, 12),1],
        NA,
        dns_model[[1]]$fcast[idxx + 5, c(3, 6, 9, 12),1],
        NA,
        dns_model[[1]]$fcast[idxx + 6, c(3, 6, 9, 12), 1]
        )

# yt <- c(dns_model[[1]]$fcast_tilt[idxx + 2, c(3, 6, 9, 12),1],
#         NA,
#         dns_model[[1]]$fcast_tilt[idxx + 3, c(3, 6, 9, 12),1],
#         NA,
#         dns_model[[1]]$fcast_tilt[idxx + 4, c(3, 6, 9, 12),1],
#         NA,
#         dns_model[[1]]$fcast_tilt[idxx + 5, c(3, 6, 9, 12),1],
#         NA,
#         dns_model[[1]]$fcast_tilt[idxx + 6, c(3, 6, 9, 12),1],
#         NA,
#         dns_model[[1]]$fcast_tilt[idxx + 7, c(3, 6, 9, 12), 1]
#         )


dd <- c(as.numeric(yield_data[seq(321,330,3),1]),
        NA,
        as.numeric(yield_data[seq(322,331,3),1]),
        NA,
        as.numeric(yield_data[seq(323,332,3),1]),
        NA,
        as.numeric(yield_data[seq(324,333,3),1]),
        NA,
        as.numeric(yield_data[seq(325,334,3),1]),
        NA,
        as.numeric(yield_data[seq(326,335,3),1])
        )

ee <- c(unclass(bc_3m_release[idxx + 180]),
        NA,
        unclass(bc_3m_release[idxx + 181]),
        NA,
        unclass(bc_3m_release[idxx + 182]),
        NA,
        unclass(bc_3m_release[idxx + 183]),
        NA,
        unclass(bc_3m_release[idxx + 184]),
        NA,
        unclass(bc_3m_release[idxx + 185]))


# ee <- c(
#   as.numeric(
#     c(
#       bc_3m_target[[1]][idxx + 180],
#       bc_3m_target[[2]][idxx + 180],
#       bc_3m_target[[3]][idxx + 180],
#       bc_3m_target[[4]][idxx + 180]
#     )
#   ),
#   NA,
#   as.numeric(
#     c(
#       bc_3m_target[[1]][idxx + 181],
#       bc_3m_target[[2]][idxx + 181],
#       bc_3m_target[[3]][idxx + 181],
#       bc_3m_target[[4]][idxx + 181]
#     )
#   ),
#   NA,
#   as.numeric(
#     c(
#       bc_3m_target[[1]][idxx + 182],
#       bc_3m_target[[2]][idxx + 182],
#       bc_3m_target[[3]][idxx + 182],
#       bc_3m_target[[4]][idxx + 182]
#     )
#   ),
#   NA,
#   as.numeric(
#     c(
#       bc_3m_target[[1]][idxx + 183],
#       bc_3m_target[[2]][idxx + 183],
#       bc_3m_target[[3]][idxx + 183],
#       bc_3m_target[[4]][idxx + 183]
#     )
#   ),
#   NA,
#   as.numeric(
#     c(
#       bc_3m_target[[1]][idxx + 184],
#       bc_3m_target[[2]][idxx + 184],
#       bc_3m_target[[3]][idxx + 184],
#       bc_3m_target[[4]][idxx + 184]
#     )
#   ),
#   NA,
#   as.numeric(
#     c(
#       bc_3m_target[[1]][idxx + 185],
#       bc_3m_target[[2]][idxx + 185],
#       bc_3m_target[[3]][idxx + 185],
#       bc_3m_target[[4]][idxx + 185]
#     )
#   )
# )

plot(
  x = xc,
  y = ee,
  type = "l",
  ylim = c(0, 1.2),
  col = "blue",
  ylab = "Yield",
  xlab = "",
  axes = F)

points(x = xc, y = ee, pch = 15, col = "blue")

lines(x = xc, y = yy, col = 'red')
points(x = xc, y = yy, col = 'red', pch = 13)

# lines(x = xc, y = dd, col = 'black')
# points(x = xc, y = dd, col = 'black', pch = 13)

# lines(x = xc, y = yt, col = 'darkgreen')
# points(x = xc, y = yt, col = 'darkgreen', pch = 13)


axis(1, labels = rep(c("1", "2", "3", "4"), 6),
     at = c(1:4, 6:9, 11:14, 16:19, 21:24, 26:29),
     cex.axis = 0.6)


mtext("June 2011", at = c(2.5), side = 1, line = 2, cex = 0.6)
mtext("July 2011", at = c(7.5), side = 1, line = 2, cex = 0.6)
mtext("August 2011", at = c(12.5), side = 1, line = 2, cex = 0.6)
mtext("September 2011", at = c(17.5), side = 1, line = 2, cex = 0.6)
mtext("October 2011", at = c(22.5), side = 1, line = 2, cex = 0.6)
mtext("November 2011", at = c(27.5), side = 1, line = 2, cex = 0.6)

axis(2)
box()
grid()
legend(x = 23, y = 1.2, c("BC Survey", "DNS Model"), col = c("blue", "red"), pch = c(15, 13), lty = c(1, 1), bty = "n", cex = 0.8,  y.intersp = 1.5)



```



```{r fomc_tilted, fig.cap="The figure shows the 12-month-ahead yield curve forecast implied by the DNS model and the corresponding anchored forecast made before and after the FOMC State- ment of August 9, 2011, together with the actual yield curve realization."}

t <- 140  ## August
m <- c(3,6,12,24,36,60)
h <- 12
idx <- match(m, maturities)

par(
  mar = c(3, 3, 2, 1),
  mgp = c(2, .7, 0),
  tck = -.01,
  cex.lab = 0.75,
  cex.axis = 0.75
)

ml <- matrix(c(1,2,3,3), nrow = 2, ncol = 2, byrow = TRUE)
layout(mat = ml, heights = c(0.45,0.1))


plot(
  as.numeric(yield_data[t + 180 + h - 1, idx]),
  lwd = 2,
  ylim = c(0, 1.5),
  lty = 3,
  pch = 1,
  type = "b",
  axes = F,
  ylab = "Yield",
  xlab = "Maturity (months)",
  main  = "August 2011",
  cex.main = 0.7
  )

axis(side = 1, labels = m, at = 1:length(m))
box()
axis(side = 2)
lines(dns_model[[1]]$fcast[t,h,idx], lty = 3, col = "red", lwd = 1.5)
points(dns_model[[1]]$fcast[t,h,idx], pch = 3, col = "red")

lines(dns_model[[1]]$fcast_tilt[t, h, idx], col = 'blue', lty = 3, lwd = 1.5)
points(dns_model[[1]]$fcast_tilt[t,h,idx], pch = 4, col = "blue")
grid()

t <- 141  ## September
plot(
  as.numeric(yield_data[t + 180 + h - 1, idx]),
  lwd = 2,
  ylim = c(0, 1.5),
  lty = 3,
  pch = 1,
  type = "b",
  ylab = "Yield",
  xlab = "Maturity (months)",
  main  = "September 2011",
  cex.main = 0.7,
  axes = F
  )
  axis(side = 1, labels = m, at = 1:length(m))
box()
axis(side = 2)


lines(dns_model[[1]]$fcast[t,h,idx], lty = 3, col = "red", lwd = 1.5)
points(dns_model[[1]]$fcast[t,h,idx], pch = 3, col = "red")

lines(dns_model[[1]]$fcast_tilt[t, h, idx], col = 'blue', lty = 3, lwd = 1.5)
points(dns_model[[1]]$fcast_tilt[t, h, idx], pch = 4, col = "blue")
grid()


# mtext("Yield",
#       side = 2, padj = 2, outer = TRUE, cex = 0.75, adj = 1)
# mtext("Maturities (months)", side = 1, padj = -2, outer = TRUE, cex = 0.75)

par(mar = c(0, 0, 1, 0))
plot(
  1,
  type = "n",
  axes = FALSE,
  xlab = "",
  ylab = "",
  frame.plot = FALSE
)
plot_colors <- c("black", "red", "blue")
legend(x = "top", inset = 0,
       legend = c("Data", "DNS", "Anchored DNS"),
       col = plot_colors, lty = c(3,3,3), pch = c(1,3,4),
       lwd = 1.5, cex = .8, horiz = TRUE, bty = "n", seg.len = 4)

```


```{r 3dyield, dev='pdf', dev.args=list(bg='transparent'), tidy=FALSE, message=FALSE, warning=FALSE, echo=FALSE, fig.width = 7, fig.height=7.5}
idix <- axTicksByTime(yield_data)
idix <- idix[seq(1, length(idix), 4)]
labels <- c("1985", "1989", "1993", "1997", "2001", "2005", "2009", "2012")
wireframe(as.matrix(yield_data),
          aspect = c(1.1, .99),
          screen = list(x=98, z = 190, y = 130),
          xlab = 'Time', ylab = 'Maturity (months)', zlab = 'Yields',
          shade = FALSE,
          scales = list(arrows = FALSE, col = 'black',
                        zlim = c(0,12.4),  
                        x = list(cex = 0.55, at = idix,
                                 labels = names(idix)),
                        y = list(cex = 0.5),
                        z = list(cex = 0.6),
                        axis=list(text=list(cex=.7))),
          par.settings = list(axis.line = list(col = "transparent"),
                              layout.heights=list(top.padding=-1, bottom.padding=-2)))
```


\clearpage

## Replication information

This file is produced from the `R` markdown document `kalman.Rmd`.

This pdf was produced on
```{r}
timestamp()
```

The specs for the platform and the `R` installation were:

```{r, echo=TRUE}
R.version
```

The following packages were used:

```{r}
for (j in pkg_from_cran) {cat(j); print(packageVersion(j))}
```
