---
title: "Average relative accuracy indices"
output: rmarkdown::html_vignette
author: Daniele Girolimetto
date: "`r Sys.Date()`"
vignette: >
  %\VignetteIndexEntry{Average relative accuracy indices}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

Let
$$
\hat{e}_{i,j,t}^{[k],h} = y_{i,t+h}^{[k]} - \hat{y}_{i,j,t}^{[k],h}, \quad
\begin{array}{l}
i=1,\ldots,n, \\
j=0,\ldots,J,
\end{array} \;  t=1,\ldots, q , \;
\begin{array}{l}
k \in {\cal K}, \\
h=1,\ldots,h_k ,
\end{array}
$$
be the forecast error, where $y$ and $\hat{y}$ are the actual and the forecasted values, respectively, suffix $i$ denotes the variable of interest, $j$ is the forecasting technique, where $j=0$ is the benchmark forecasting procedure, $t$ is the forecast origin, ${\cal K}$ is the set of the time frequencies at which the series is observed, and $h$ is the forecast horizon, whose lead time depends on the time frequency $k$.

Denote by A$_{i,j}^{[k],h}$ the forecasting accuracy of the technique $j$, computed across $q$ forecast origins, for the $h$-step-ahead forecasts of the variable $i$ at the temporal aggregation level $k$. For example, A$_{i,j}^{[k],h} = MSE_{i,j}^{[k],h}$ (*type="mse"* in `score_index()` function, *default*), otherwise we might have A$_{i,j}^{[k],h} = MAE_{i,j}^{[k],h}$ (*type="mae"* in `score_index()` function) or A$_{i,j}^{[k],h} = RMSE_{i,j}^{[k],h}$ (*type="rmse"* in `score_index()` function), where
$$
\begin{array}{rcl}
MSE_{i,j}^{[k],h} & = & \displaystyle\frac{1}{q}\displaystyle\sum_{t=1}^{q} \left(\hat{e}_{i,j}^{[k],h}\right)^2 \\
MAE_{i,j}^{[k],h} & = & \displaystyle\frac{1}{q}\displaystyle\sum_{t=1}^{q} \left|\hat{e}_{i,j}^{[k],h}\right| \\
RMSE_{i,j}^{[k],h} & = & \sqrt{\displaystyle\frac{1}{q}\displaystyle\sum_{t=1}^{q} \left(\hat{e}_{i,j}^{[k],h}\right)^2}
\end{array}
$$

\noindent In any case, we consider the relative version of the accuracy index $A_{i,j}^{[k],h}$, given by:
$$
r_{i,j}^{[k],h} = \displaystyle\frac{A_{i,j}^{[k],h}}{A_{i,0}^{[k],h}} ,
\quad i=1,\ldots,n, \quad  j=0,\ldots, J, \quad k \in {\cal K}, \quad h=1, \ldots, h_k,
$$
and use it to compute the Average relative accuracy index of the forecasting procedure $j$, for given $k$ and $h$, through the geometric mean:
$$
\text{AvgRelA}_{j}^{[k],h} = \left(\displaystyle\prod_{i=1}^{n} r_{i,j}^{[k],h} \right)^{\frac{1}{n}} , \quad j=0,\ldots,J .
$$

\noindent We may consider the following average relative accuracy indices for selected groups of variables/time frequencies and forecast horizons:

**Average relative accuracy indices for a single variable at a given time frequency, for multiple forecast horizons**
$$
\label{AvgRelAikh1h2}
\text{AvgRelA}_{i,j}^{[k],q_1:q_2} = \left(\prod_{h=q_1}^{q_2} r_{i,j}^{[k],h}\right)^{\frac{1}{q_2 - q_1 + 1}}, \;
\begin{array}{l}
	i=1,\ldots,n, \\
	j=0,\ldots,J,
\end{array} \;
\begin{array}{l}
	k \in {\cal K}, \\
	1 \le q_1 \le q_2 \le h_k
\end{array} .
%i=1,\ldots,n , \; j=0,\ldots, J, \; k \in {\cal K}, \; 1 \le q_1 \le q_2 \le h_k.
$$
**Average relative accuracy indices for a group of variables (either all, or selected groups, e.g. a: uts, b: bts) at a given time frequency, either for a single forecast horizon or across them**
$$
\begin{array}{rcll}
\text{AvgRelA}^{[k],h}_j & = &
\left(\displaystyle\prod_{i=1}^{n} r_{i,j}^{[k],h}\right)^{\frac{1}{n}}, &
 j=0,\ldots,J , \; k \in {\cal K}, \; h=1,\ldots,h_k \\[.5cm]
\text{AvgRelA}^{[k],h}_{a,j} & = & \left(\displaystyle\prod_{i=1}^{n_a} r_{i,j}^{[k],h}\right)^{\frac{1}{n_a}}, &
 j=0,\ldots, J, \; k \in {\cal K} \\[.15cm]
\text{AvgRelA}^{[k],h}_{b,j} & = & \left(\displaystyle\prod_{i=n_a+1}^{n} r_{i,j}^{[k],h}\right)^{\frac{1}{n_b}}, &
j=0,\ldots, J, \; k \in {\cal K} \\[.15cm]
\text{AvgRelA}^{[k]}_j & = & \left(\displaystyle\prod_{i=1}^{n} \prod_{h=1}^{h_k} r_{i,j}^{[k],h}\right)^{\frac{1}{n h_k}}, &
\; j=0,\ldots,J, \; k \in {\cal K} \\[.15cm]
\text{AvgRelA}^{[k]}_{a,j} & = &\left(\displaystyle\prod_{i=1}^{n_a} \prod_{h=1}^{h_k} r_{i,j}^{[k],h}\right)^{\frac{1}{n_a h_k}}, &
j=0,\ldots, J, \; k \in {\cal K} \\[.15cm]
\text{AvgRelA}^{[k]}_{b,j} & = &\left(\displaystyle\prod_{i=n_a+1}^{n} \prod_{h=1}^{h_k} r_{i,j}^{[k],h}\right)^{\frac{1}{n_b h_k}}, &
 j=0,\ldots, J, \; k \in {\cal K}
\end{array}
$$

**Average relative accuracy indices for a single variable or for a group of variables (all, a: uts, b: bts), across all time frequencies and forecast horizons**
$$
\begin{array}{rcll}
\text{AvgRelA}_{i,j} & = & \left(\displaystyle\prod_{k \in {\cal K}} \prod_{h=1}^{h_k} r_{i,j}^{[k],h}\right)^{\frac{1}{k^*+m}}, &
\begin{array}{l}
i=1,\ldots,n  \\ j=0,\ldots,J
\end{array} \\[.15cm]
\text{AvgRelA}_j & = & \left(\displaystyle\prod_{i=1}^{n} \prod_{k \in {\cal K}} \prod_{h=1}^{h_k} r_{i,j}^{[k],h}\right)^{\frac{1}{n(k^*+m)}}  , & j=0,\ldots,J \\[.15cm]
\text{AvgRelA}_{a,j} & = & \left(\displaystyle\prod_{i=1}^{n_a} \prod_{k \in {\cal K}} \prod_{h=1}^{h_k} r_{i,j}^{[k],h}\right)^{\frac{1}{n_a(k^*+m)}}  , & j=0,\ldots,J \\[.15cm]
\text{AvgRelA}_{b,j} & = & \left(\displaystyle\prod_{i=n_a+1}^{n} \prod_{k \in {\cal K}} \prod_{h=1}^{h_k} r_{i,j}^{[k],h}\right)^{\frac{1}{n_b(k^*+m)}} , & j=0,\ldots,J
\end{array}
$$

## The `score_index()` function outputs

The `score_index()` function returns a summary table called **Avg$\_$mat** (if *compact* option is **TRUE**, *default*), otherwise it returns a list of four tables. Fixed method $j$, we have:

- **Avg$\_$mat** (*default* output)
```{r table-simple, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
df <- cbind(all = c("$\\textbf{all}$","$\\text{AvgRelA}^{[m]}_j$","$\\vdots$","$\\text{AvgRelA}^{[k]}_j$","$\\vdots$",
                         "$\\text{AvgRelA}^{[1]}_j$","$\\vdots$","$\\text{AvgRelA}_j$"),
                 uts = c("$\\textbf{uts}$","$\\text{AvgRelA}^{[m]}_{a,j}$","$\\vdots$","$\\text{AvgRelA}^{[k]}_{a,j}$","$\\vdots$",
                         "$\\text{AvgRelA}^{[1]}_{a,j}$","$\\vdots$","$\\text{AvgRelA}_{a,j}$"),
                 bts = c("$\\textbf{bts}$","$\\text{AvgRelA}^{[m]}_{b,j}$","$\\vdots$","$\\text{AvgRelA}^{[k]}_{b,j}$","$\\vdots$",
                         "$\\text{AvgRelA}^{[1]}_{b,j}$","$\\vdots$","$\\text{AvgRelA}_{b,j}$"))
rownames(df) <- c("","$\\textbf{m}$","$\\vdots$","$\\textbf{k}$","$\\vdots$","$\\textbf{1}$","$\\vdots$","$\\textbf{all}$")
knitr::kable(df,align='cccc',escape = F, col.names = rep("",3))
```

- **Rel$\_$mat** (if *compact = FALSE*)
```{r table-simple2, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
df <- cbind(col1 = c("$\\textbf{m}$","$\\textbf{1}$","$\\text{RelA}^{[m],1}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[m],1}_{i,j}$","$\\vdots$","$\\text{RelA}^{[m],1}_{n,j}$"),
            col2 = c("$\\dots$","$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col3 = c("", "$\\textbf{1}$","$\\text{RelA}^{[k],1}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[k],1}_{i,j}$","$\\vdots$","$\\text{RelA}^{[k],1}_{n,j}$"),
            col4 = c("$\\textbf{k}$","$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col5 = c("", "$\\mathbf{h_k}$","$\\text{RelA}^{[k],h_k}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[k],h_k}_{i,j}$","$\\vdots$","$\\text{RelA}^{[k],h_k}_{n,j}$"),
            col6 = c("$\\dots$","$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col7 = c("", "$\\textbf{1}$","$\\text{RelA}^{[1],1}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[1],1}_{i,j}$","$\\vdots$","$\\text{RelA}^{[1],1}_{n,j}$"),
            col8 = c("$\\textbf{1}$","$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col9 = c("", "$\\mathbf{m}$","$\\text{RelA}^{[1],m}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[1],m}_{i,j}$","$\\vdots$","$\\text{RelA}^{[1],m}_{n,j}$"))
rownames(df) <- c("${\\cal K}$","$\\textbf{h}$","$\\textbf{1}$","$\\vdots$","$\\textbf{i}$","$\\vdots$","$\\textbf{n}$")
knitr::kable(df,align='ccccccccc',escape = F, col.names = rep("",9))
```

- **Rel$\_$mat$\_$cum ** (if *compact = FALSE*)
```{r table-simple2cum, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
df <- cbind(col1 = c("$\\textbf{m}$","$\\textbf{1:1}$","$\\text{RelA}^{[m],1:1}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[m],1:1}_{i,j}$","$\\vdots$","$\\text{RelA}^{[m],1:1}_{n,j}$"),
            col2 = c("$\\dots$","$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col3 = c("", "$\\textbf{1:1}$","$\\text{RelA}^{[k],1:1}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[k],1:1}_{i,j}$","$\\vdots$","$\\text{RelA}^{[k],1:1}_{n,j}$"),
            col4 = c("$\\textbf{k}$","$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col5 = c("", "$\\mathbf{1:h_k}$","$\\text{RelA}^{[k],1:h_k}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[k],1:h_k}_{i,j}$","$\\vdots$","$\\text{RelA}^{[k],1:h_k}_{n,j}$"),
            col6 = c("$\\dots$","$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col7 = c("", "$\\textbf{1:1}$","$\\text{RelA}^{[1],1:1}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[1],1:1}_{i,j}$","$\\vdots$","$\\text{RelA}^{[1],1:1}_{n,j}$"),
            col8 = c("$\\textbf{1}$","$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col9 = c("", "$\\mathbf{1:m}$","$\\text{RelA}^{[1],1:m}_{1,j}$","$\\vdots$",
                         "$\\text{RelA}^{[1],1:m}_{i,j}$","$\\vdots$","$\\text{RelA}^{[1],1:m}_{n,j}$"))
rownames(df) <- c("${\\cal K}$","$\\textbf{h}$","$\\textbf{1}$","$\\vdots$","$\\textbf{i}$","$\\vdots$","$\\textbf{n}$")
knitr::kable(df,align='ccccccccc',escape = F, col.names = rep("",9))
```
- **AvgRelIND$\_$ik** (if *compact = FALSE*)
```{r table-simple3, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
df <- cbind(col1 = c("$\\textbf{m}$","$\\text{AvgRelA}^{[m]}_{1,j}$","$\\vdots$",
                         "$\\text{AvgRelA}^{[m]}_{i,j}$","$\\vdots$","$\\text{AvgRelA}^{[m]}_{n,j}$"),
            col2 = c("$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col3 = c("$\\textbf{k}$", "$\\text{AvgRelA}^{[k]}_{1,j}$","$\\vdots$",
                         "$\\text{AvgRelA}^{[k]}_{i,j}$","$\\vdots$","$\\text{AvgRelA}^{[k]}_{n,j}$"),
            col4 = c("$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col5 = c("$\\textbf{1}$","$\\text{AvgRelA}^{[1]}_{1,j}$","$\\vdots$",
                         "$\\text{AvgRelA}^{[1]}_{i,j}$","$\\vdots$","$\\text{AvgRelA}^{[1]}_{n,j}$"),
            col6 = c("$\\dots$","$\\dots$","","$\\dots$","","$\\dots$"),
            col7 = c("$\\mathbf{all}$","$\\text{AvgRelA}_{1,j}$","$\\vdots$",
                         "$\\text{AvgRelA}_{i,j}$","$\\vdots$","$\\text{AvgRelA}_{n,j}$"))
rownames(df) <- c("","$\\textbf{1}$","$\\vdots$","$\\textbf{i}$","$\\vdots$","$\\textbf{n}$")
knitr::kable(df,align='ccccccc',escape = F, col.names = rep("",7))
```	
	
- **Avg$\_$k** (if *compact = FALSE*)
```{r table-simple4, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
df <- cbind(col1 = c("$\\textbf{m}$","$\\textbf{1}$","$\\text{AvgRelA}^{[m],1}_{j}$",
                         "$\\text{AvgRelA}^{[m],1}_{a,j}$","$\\text{AvgRelA}^{[m],1}_{b,j}$"),
            col2 = c("$\\dots$","$\\dots$","$\\dots$","$\\dots$","$\\dots$"),
            col3 = c("", "$\\textbf{1}$","$\\text{AvgRelA}^{[k],1}_{j}$",
                         "$\\text{AvgRelA}^{[k],1}_{a,j}$","$\\text{AvgRelA}^{[k],1}_{b,j}$"),
            col4 = c("$\\textbf{k}$","$\\dots$","$\\dots$","$\\dots$","$\\dots$"),
            col5 = c("", "$\\mathbf{h_k}$","$\\text{AvgRelA}^{[k],h_k}_{j}$",
                         "$\\text{AvgRelA}^{[k],h_k}_{a,j}$","$\\text{AvgRelA}^{[k],h_k}_{b,j}$"),
            col6 = c("$\\dots$","$\\dots$","$\\dots$","$\\dots$","$\\dots$"),
            col7 = c("", "$\\textbf{1}$","$\\text{AvgRelA}^{[1],1}_{j}$",
                         "$\\text{AvgRelA}^{[1],1}_{a,j}$","$\\text{AvgRelA}^{[1],1}_{b,j}$"),
            col8 = c("$\\textbf{1}$","$\\dots$","$\\dots$","$\\dots$","$\\dots$"),
            col9 = c("", "$\\mathbf{m}$","$\\text{AvgRelA}^{[1],m}_{j}$",
                         "$\\text{AvgRelA}^{[1],m}_{a,j}$","$\\text{AvgRelA}^{[1],m}_{b,j}$"))
rownames(df) <- c("${\\cal K}$","$\\textbf{h}$","$\\textbf{all}$","$\\textbf{a}$","$\\textbf{b}$")
knitr::kable(df,align='ccccccccc',escape = F, col.names = rep("",9))
```

- **Avg$\_$k$\_$cum** (if *compact = FALSE*)
```{r table-simple4cum, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
df <- cbind(col1 = c("$\\textbf{m}$","$\\textbf{1:1}$","$\\text{AvgRelA}^{[m],1:1}_{j}$",
                         "$\\text{AvgRelA}^{[m],1:1}_{a,j}$","$\\text{AvgRelA}^{[m],1:1}_{b,j}$"),
            col2 = c("$\\dots$","$\\dots$","$\\dots$","$\\dots$","$\\dots$"),
            col3 = c("", "$\\textbf{1:1}$","$\\text{AvgRelA}^{[k],1:1}_{j}$",
                         "$\\text{AvgRelA}^{[k],1:1}_{a,j}$","$\\text{AvgRelA}^{[k],1:1}_{b,j}$"),
            col4 = c("$\\textbf{k}$","$\\dots$","$\\dots$","$\\dots$","$\\dots$"),
            col5 = c("", "$\\mathbf{1:h_k}$","$\\text{AvgRelA}^{[k],1:h_k}_{j}$",
                         "$\\text{AvgRelA}^{[k],1:h_k}_{a,j}$","$\\text{AvgRelA}^{[k],1:h_k}_{b,j}$"),
            col6 = c("$\\dots$","$\\dots$","$\\dots$","$\\dots$","$\\dots$"),
            col7 = c("", "$\\textbf{1:1}$","$\\text{AvgRelA}^{[1],1:1}_{j}$",
                         "$\\text{AvgRelA}^{[1],1:1}_{a,j}$","$\\text{AvgRelA}^{[1],1:1}_{b,j}$"),
            col8 = c("$\\textbf{1}$","$\\dots$","$\\dots$","$\\dots$","$\\dots$"),
            col9 = c("", "$\\mathbf{1:m}$","$\\text{AvgRelA}^{[1],1:m}_{j}$",
                         "$\\text{AvgRelA}^{[1],1:m}_{a,j}$","$\\text{AvgRelA}^{[1],1:m}_{b,j}$"))
rownames(df) <- c("${\\cal K}$","$\\textbf{h}$","$\\textbf{all}$","$\\textbf{a}$","$\\textbf{b}$")
knitr::kable(df,align='ccccccccc',escape = F, col.names = rep("",9))
```

### Examples

```{r eval=FALSE}
library(FoReco)
data(FoReco_data)

# Cross-temporal framework
oct_recf <- octrec(FoReco_data$base, m = 12, C = FoReco_data$C,
                   comb = "bdshr", res = FoReco_data$res)$recf
oct_score <- score_index(recf = oct_recf,
                         base = FoReco_data$base,
                         test = FoReco_data$test, m = 12, nb = 5)

# Cross-sectional framework#'
# monthly base forecasts
id <- which(simplify2array(strsplit(colnames(FoReco_data$base), split = "_"))[1, ] == "k1")
mbase <- t(FoReco_data$base[, id])
# monthly test set
mtest <- t(FoReco_data$test[, id])
# monthly residuals
id <- which(simplify2array(strsplit(colnames(FoReco_data$res), split = "_"))[1, ] == "k1")
mres <- t(FoReco_data$res[, id])
# monthly reconciled forecasts
mrecf <- htsrec(mbase, C = FoReco_data$C, comb = "shr", res = mres)$recf
# score
hts_score <- score_index(recf = mrecf, base = mbase, test = mtest, nb = 5)

# Temporal framework
data(FoReco_data)
# top ts base forecasts ([lowest_freq' ...  highest_freq']')
topbase <- FoReco_data$base[1, ]
# top ts residuals ([lowest_freq' ...  highest_freq']')
topres <- FoReco_data$res[1, ]
# top ts test ([lowest_freq' ...  highest_freq']')
toptest <- FoReco_data$test[1, ]
# top ts recf ([lowest_freq' ...  highest_freq']')
toprecf <- thfrec(topbase, m = 12, comb = "acov", res = topres)$recf
# score
thf_score <- score_index(recf = toprecf, base = topbase, test = toptest, m = 12)
```

## References
Di Fonzo, T., Girolimetto, D. (2021), Cross-temporal forecast reconciliation: Optimal combination method and heuristic alternatives, _International Journal of Forecasting_, in press.
