dissertacao.qmd

---
format:
  pdf:
    include-in-header: [config/tema/preamble.tex, config/tema/customizacao.tex]
    include-before-body: config/elementos/pre_textuais.tex
    include-after-body: config/elementos/pos_textuais.tex
    documentclass: abntex2
    classoption: [12pt, oneside, a4paper, chapter=TITLE, section=TITLE, brazil]
    keep-tex: false
    latex-max-runs: 4
    output-file: dissertacao.pdf
---

```{r config}
#| include = FALSE

# opções
knitr::opts_chunk$set(
  out.width = "90%"
)

# reprodutibilidade
set.seed(1)

# pacotes
pacman::p_load(
  kableExtra,
  ggplot2,
  tsibble,
  fable,
  fabletools
)

# tema ggplot
tema = theme(
  text = element_text(family = "serif"),
  legend.background = element_rect(fill = "transparent")
  )

# gerar bibliografia de pacotes
knitr::write_bib(
  c(
    .packages(),
    "basedosdados",
    "mlr3",
    "mlr3mbo",
    "xgboost",
    "ranger",
    "lightgbm",
    "glmnet",
    "hts"
  ),
  file = "config/elementos/packages.bib"
)

# dados
estban = readRDS("data/estban/estban.rds")

# filtrando dados
estban = subset(estban,
  select = c(
    ref,
    nome_mesorregiao,
    nome_microrregiao,
    nome,
    cnpj_agencia,
    verbete,
    saldo
  )
) |>
  tsibble::filter_index(~"2022-12-31") |>
  (\(x) x[order(
    x$ref,
    x$nome_mesorregiao,
    x$nome_microrregiao,
    x$nome,
    x$cnpj_agencia,
    x$verbete
  ), ])()
```

# INTRODUÇÃO

## Contextualização da pesquisa

Diversas indústrias se enquadram em categorias de negócio que requerem previsões de múltiplas séries temporais correlacionadas que são resultados de agregação. Por exemplo, o total de empréstimos de uma instituição financeira corresponde ao agregado dos empréstimos de cada uma de suas agências; o total de vendas de uma rede nacional de farmácias corresponde ao agregado de vendas de suas unidades em cada estado; o total da produção de uma petrolífica multinacional corresponde ao total produzido em cada país por cada uma de suas plataformas. A essas estruturas naturais de agregação dá-se o nome de *séries temporais hierárquicas* [@hyndman_forecasting_2021].

No que concerne a elaboração dessas previsões, pode-se realizar previsões individualmente para todos os níveis da estrutura. No caso de uma empresa qualquer, isso significa realizar previsões, por exemplo, para cada filial, para o agregado das filiais de cada região e para o agregado da instituição. Infelizmente, não há qualquer razão matemática, exceto para métodos de previsão muito simples, para que essas previsões sejam *coerentes* (i.e. que a soma das previsões individuais seja igual às previsões do agregado). Para garantir previsões coerentes entre si foram desenvolvidos os chamados métodos de previsão hierárquicas [@athanasopoulos_hierarchical_2009], sendo os mais simples o *top-down*, *bottom-up* e uma combinação das duas, o *middle-out*.

Ao menos desde a década de 1960, diversos trabalhos foram publicadas acerca das abordagens *bottom-up* e *top-down*, suas vantagens e desvantagens, e tentativas de se definir qual é o melhor método^[Revisões sistemáticas das literaturas iniciais relacionadas ao tema podem ser encontradas em @athanasopoulos_hierarchical_2009 e @athanasopoulos_forecast_2023.]. A abordagem *top-down* envolve realizar previsões para o total e então distribuí-las para cada nó filho da estrutura seguindo alguma lógica proporcional. Tomando o caminho inverso, a abordagem *bottom-up* consiste em realizar previsões para cada série do nível mais desagregado e, então, agregá-las para obter as previsões para os níveis superiores. Ambos métodos possuem suas vantagens e desvantagens e cabe ao analista avaliar o *trade-off* entre os ganhos de precisão percebidos com a geração de previsões individuais para todas as séries do menor nível hierárquico e a economia de tempo e processamento em realizar o contrário [@gross_disaggregation_1990].

Além disso, ambas são abordagens de nível único, isto é, são realizadas as previsões para um único nível e então os demais níveis são obtidos agregando ou desagregando. O problema com esses tipos de abordagem é que elas utilizam informação incompleta [@hyndman_forecasting_2021]. Por exemplo, suponha-se que se escolha estimar modelos para cada uma das $n$ filiais de uma empresa e agregá-las (*bottom-up*). Nesse caso, ignora-se a influência que os níveis mais agregados podem ter na estimação de cada filial. Por outro lado, caso se escolha estimar modelos para o nível mais agregado (*top-down*), ignora-se a informação individual de cada filial.

Uma terceira possibilidade é a *reconciliação ótima*. Ela é uma abordagem que busca resolver esse problema e consiste em realizar previsões para todos os níveis hierárquicos e, então, estimar um modelo para reescrever as previsões do nível mais desagregado como uma combinação linear de todos os elementos da hierarquia. Obtidas as novas previsões no menor nível, elas são então agregadas, gerando previsões coerentes nos níveis superiores. Dessa forma, a informação de todos os níveis é utilizada na estimação dos modelos e na geração das previsões ao mesmo tempo em que a variância do erro de previsão é minimizada [@wickramasuriya_optimal_2019].

Atualmente, os métodos analíticos baseados na minimização do traço (MinT) da matriz da variância-covariância dos erros das previsões reconciliadas, desenvolvidos em  @wickramasuriya_optimal_2019, estão entre os mais populares na literatura da reconciliação ótima. Esses métodos divergem apenas na forma da qual se dá o relacionamento entre os diferentes elementos da hierarquia, resultando numa estimação por mínimos quadrados ordinários (MQO), mínimos quadrados ponderados (MQP) ou mínimos quadrados generalizados (MQG)^[Se os erros de previsão são descorrelacionados e homoscedásticos ao longo de toda estrutura (MQO), o que é impossível em séries temporais hierárquicas; se os erros são descorrelacionados e homoscedásticos apenas dentro do mesmo nível hierárquico (MQP estrutural); se os erros são descorrelacionados mas ponderados pela variância da série (MQP); ou se são correlacionados e variantes ao longo de toda a estrutura (MQG).].

Entretanto, tais métodos são sujeitos a uma série de restrições, como as do modelo clássico de regressão linear (MCRL), e têm sua capacidade preditiva reduzida quando suas hipóteses são violadas. Além disso, esses métodos requerem que *toda* a informação seja utilizada, mesmo aquelas que eventualmente não sejam relevantes para aquela previsão. Nesse sentido, métodos de *machine learning* são mais gerais, de modo a permitir parâmetros não lineares e poderem aproximar virtualmente qualquer função, além de incluírem parâmetros de regularização, de forma que não necessariamente utilizam toda a informação disponível. Espera-se, portanto, que esses métodos alcancem melhor performance no problema da reconciliação ótima, justificando a pesquisa e atenção ao tema. Nessa direção, @spiliotis_hierarchical_2021 propuseram uma metodologia para reconciliação ótima de séries temporais hierárquicas utilizando métodos de *machine learning*, especificamente o *gradient boosting* e *random forest*, obtendo resultados superiores aos métodos analíticos em algumas situações.

Tomando como ponto de partida a metodologia proposta por @spiliotis_hierarchical_2021, este trabalho explora métodos adicionais de ML na tarefa de reconciliação ótima, especificamente os métodos de regressão regularizada Lasso, *ridge* e *elastic net*, e *support vector machines* (SVM), além de outra implementação de *gradient boosting*, o *LightGBM* (LGBM). Paralelamente, propõe e avalia estratégias alternativas na obtenção do conjunto de treinamento para os modelos de ML, verificando se sua performance preditiva supera a dos métodos analíticos nas situações em que a metodologia original não foi capaz.

Para tanto, foram realizados três estudos de caso, sendo o primeiro em um contexto de economia bancária brasileira, enquanto os demais envolvem bases de dados amplamente utilizadas na literatura de séries temporais hierárquicas, relacionadas ao setor de turismo doméstico australiano. O primeiro estudo de caso trata da previsão de saldos de empréstimos e financiamentos do Banco do Estado do Espírito Santo (Banestes); o segundo, da previsão da quantidade trimestral de viagens domésticas na Austrália; e o terceiro, da previsão da quantidade mensal de noites que os visitantes passam em viagens domésticas na Austrália.

## Outline

Esta dissertação está organizada da seguinte forma: no Capítulo 2, é introduzida a formalização algébrica e metodológica necessária para a compreensão dos métodos de previsões de séries temporais hierárquicas e reconciliação ótima, assim como uma breve revisão de literatura englobando os artigos seminais para o tema. No Capítulo 3, são apresentadas as propostas de variações metodológicas para combinação ótima por ML. No Capítulo 4, é realizado o estudo de caso para a previsão de saldos de crédito e financiamento para o Banestes. No Capítulo 5, são realizados os estudos de casos para as bases de turismo doméstico australiano.

Uma síntese das contribuições desse trabalho é expressa pelos seguintes pontos:

1. Como trabalhos anteriores já haviam evidenciado, a reconciliação ótima por ML é capaz de melhorar a acurácia das previsões de séries temporais hierárquicas e agrupadas em relação aos métodos analíticos;
1. Além do *XGBoost* (XGB) e do *random forest* (RF), já testados em @spiliotis_hierarchical_2021, foram avaliados na tarefa de reconciliação ótima os métodos de ML LGBM, *Lasso*, *ridge*, *elastic net*, e SVM.
1. Principalmente devido à quantidade de hiperparâmetros a serem calibrados, os métodos de ML baseados em árvore de decisão (XGB, RF e LGBM) apresentaram custo elevado em tempo de processamento, ao mesmo tempo em que não foram capazes de superar métodos mais simples, como os de regressão regularizada e SVM, de forma que o *trade-off* entre acurácia e tempo de processamento não foi justificado.
1. Foram investigadas três estratégias para obtenção do conjunto de treinamento dos modelos de reconciliação ótima por ML: o *rolling forecasting* (proposta em @spiliotis_hierarchical_2021), e outras duas adicionais propostas neste trabalho, a primeira utilizando os valores ajustados dos modelos de previsão base e a segunda utilizando os valores ajustados de modelos reestimados;
1. Para otimização dos hiperparâmetros dos modelos de ML, foi proposta a utilização de validação cruzada com origem móvel;
1. Não houve método ou estratégia que alcançasse consistentemente a melhor performance em todos os três estudos de caso, ou mesmo ao longo de todos os níveis hierárquicos de cada estudo de caso. Em geral, os métodos de ML obtiveram melhor resultados, ou nos níveis mais desagregados, ou nos mais agregados, a depender do conjunto de dados analisado;
1. A combinação correta do método de ML e da estratégia para obtenção do conjunto treinamento pode resultar em até 93% de ganho de performance em relação aos métodos analíticos de reconciliação ótima.

# RECONCILIAÇÃO DE SÉRIES TEMPORAIS HIERÁRQUICAS E AGRUPADAS

## Notação algébrica

Séries temporais hierárquicas são aquelas que podem ser agregadas ou desagregadas naturalmente em uma estrutura aninhada [@hyndman_forecasting_2021]. Para ilustrar, tome a série do PIB de um país fictício com três estados, cada um com dois municípios. Essa série pode ser desagregada por estado que, por sua vez, pode ser desagregada por município (@fig-h).

![Séries Hierárquicas](img/hierarq.png){#fig-h}

Essa estrutura pode ser representada através de equações para qualquer nível de agregação. Dessa forma, o agregado nacional pode ser descrito pelos agregados dos estados, Equação \eqref{eq:ha}, ou como o agregado dos municípios, Equação \eqref{eq:ha_mun}. Já o agregado para o estado A é representado pela Equação \eqref{eq:haES}.

\begin{align}
y_t &= y_{A,t} + y_{B,t} + y_{C,t} \label{eq:ha} \\
y_t &= y_{AA,t} + y_{AB,t} + y_{BA,t} + y_{BB,t} + y_{CA,t} + y_{CB,t}\label{eq:ha_mun} \\
y_{A,t} &= y_{AA,t} + y_{AB,t}\label{eq:haES}
\end{align}

Alternativamente, podemos descrever a estrutura completa de forma matricial:

$$
\begin{bmatrix}
    y_{t} \\
    y_{A, t} \\
    y_{B, t} \\
    y_{C, t} \\
    y_{AA, t} \\
    y_{AB, t} \\
    y_{BA, t} \\
    y_{BB, t} \\
    y_{CA, t} \\
    y_{CB, t}
\end{bmatrix}_{10 \times 1}
=
\begin{bmatrix}
    1 & 1 & 1 & 1 & 1 & 1 \\
    1 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 1 \\
    1 & 0 & 0 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}_{10 \times 6}
\begin{bmatrix}
    y_{AA, t} \\
    y_{AB, t} \\
    y_{BA, t} \\
    y_{BB, t} \\
    y_{CA, t} \\
    y_{CB, t}
\end{bmatrix}_{6 \times 1}
$$ {#eq-matriz_hierarquia}

Uma outra forma de desagregarmos o PIB é por atividade econômica --- agricultura, indústrias extrativas, indústria de transformação, eletricidade e gás, construção etc (@fig-a). Essa estrutura não pode ser desagregada naturalmente de uma única maneira, como é a hierarquia de estados e municípios. Não pode ser aninhada por um atributo como a própria geografia. A esse tipo de estrutura dá-se o nome de *séries temporais agrupadas*.

![Séries Agrupadas](img/agrupadas.png){#fig-a}

Combinando as duas, temos a estrutura de séries hierárquicas agrupadas. Ao contrário da estrutura hierárquica, que só pode ser agregada de uma forma, como com os municípios abaixo dos estados^[Essa estrutura é única no sentido que o somatório dos municípios totaliza o estado, mas não se pode somar estados para totalizar um município. Outro exemplo de estrutura hierárquica é a série de vendas de uma empresa: pode-se agregar as vendas de cada filial para obter o total, mas não o contrário.], a adição da estrutura agrupada pode ocorrer tanto acima (@fig-ha1) quanto abaixo (@fig-ha2) da hierárquica.

![Séries Hierárquicas Agrupadas (a)](img/hier_agrup.png){#fig-ha1}

![Séries Hierárquicas Agrupadas (b)](img/hier_agrup_2.png){#fig-ha2}

Na notação matricial, a estrutura da @fig-ha2 é representada como abaixo. Formalmente, o primeiro membro da igualdade é composto pelo vetor $\mathbfit{y}_t$ $n$-dimensional com todas as observações no tempo $t$ para todos os níveis da hierarquia. O segundo membro é composto pela matriz de soma $\mathbfit{S}$ de dimensão $n \times m$ que define as equações para todo nível de agregação, e pelo vetor $\mathbfit{b}_t$ composta pelas séries no nível mais desagregado.

$$
\mathbfit{y}_t=\mathbfit{Sb}_t
$$ {#eq-vetor_b}

$$
\begin{bmatrix}
    y_{t} \\
    y_{A, t} \\
    y_{B, t} \\
    y_{C, t} \\
    y_{X, t} \\
    y_{Y, t} \\
    y_{AX, t} \\
    y_{AY, t} \\
    y_{BX, t} \\
    y_{BY, t} \\
    y_{CX, t} \\
    y_{CY, t}
\end{bmatrix}_{12 \times 1}
=
\begin{bmatrix}
    1 & 1 & 1 & 1 & 1 & 1 \\
    1 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 1 \\
    1 & 0 & 1 & 0 & 1 & 0 \\
    0 & 1 & 0 & 1 & 0 & 1 \\
    1 & 0 & 0 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 0 & 0 & 1 \\
\end{bmatrix}_{12 \times 6}
\begin{bmatrix}
    y_{AX, t} \\
    y_{AY, t} \\
    y_{BX, t} \\
    y_{BY, t} \\
    y_{CX, t} \\
    y_{CY, t}
\end{bmatrix}_{6 \times 1}
$$ {#eq-matriz_ha}

## Abordagens de nível único

As abordagens de nível único são aquelas em que as previsões são realizadas para um único nível hierárquico. A partir dessas previsões, os demais níveis são obtidos, ou desagregando (no caso dos níveis inferiores), ou agregando (no caso dos níveis superiores) essas previsões [@hyndman_forecasting_2021]. Os métodos *top-down*, *bottom-up* e *middle-out* são abordagens de nível único.

Enquanto há apenas uma única forma de se agregar níveis na hierarquia (*bottom-up*), a desagregação (*top-down*) pode ser realizada de, ao menos, duas dezenas de maneiras [@gross_disaggregation_1990]. Dois dos métodos mais intuitivos são a média das proporções históricas e a proporção das médias históricas.

Na média das proporções históricas, cada proporção $p_j$, com $j = {1,...,m}$, consiste em tomar a média das proporções da série desagregada $y_{j,t}$ em relação ao agregado $y_{1,t}$:

$$
p_j = \frac{1}{T} \sum^{T}_{t=1} \frac{y_{j,t}}{y_{1,t}}
$$ {#eq-td_1}

Já a proporção das médias históricas consiste em tomar a proporção das médias das séries desagregadas em relação à média do agregado^[Isso é equivalente a tomar a proporção direta entre os somatórios das séries. Note que, pelas propriedades do operador de somatório, $\sum^{T}_{t=1} \frac{y_{t}}{T} = \frac{y_1}{T}+...+\frac{y_T}{T} = \frac{y_1+...+y_T}{T} = \frac{\sum^{T}_{t=1} y_t}{T}$. Então, a equação @eq-td_2 pode ser simplificada para $p_j = \frac{\sum^{T}_{t=1} y_{j,t}}{\sum^{T}_{t=1} y_{1,t}}$.].

$$
p_j = \frac{\sum^{T}_{t=1} \frac{y_{j,t}}{T}}{\sum^{T}_{t=1} \frac{y_{1,t}}{T}}
$$ {#eq-td_2}

@athanasopoulos_hierarchical_2009 desenvolvem o método proporções de previsão, que consiste em um método *top-down* em que os pesos são calculados a partir das proporções das previsões fora da amostra ao invés do passado. A vantagem do método é que os pesos estarão os mais próximos das características mais recentes da série, ao invés de serem baseados em dados históricos. A desvantagem é que se deve realizar previsões para toda a hierarquia, perdendo o ganho de agilidade dos demais métodos *top-down*.

$$
p_j = \prod^{K-1}_{i=0}\frac{\hat{Y}^{i}_{j,t+h}}{\sum\hat{Y}^{(i+1)}_{j,t+h}}
$$

Uma vez definida a forma e os pesos com os quais as previsões do agregado serão distribuídas nos níveis inferiores, podemos definir uma matriz com todos esses pesos, que, seguindo a formulação de @hyndman_forecasting_2021, chamamos de $\mathbfit{G}$:

$$
\mathbfit{G}
=
\begin{bmatrix}
    p_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}_{6 \times 10}
$$ {#eq-matriz_g}

$\mathbfit{G}$ é uma matriz $m \times n$ que multiplica o vetor $\hat{\mathbfit{y}}_{T+h|T}$ que chamamos de *previsões base*, isto é, as previsões individuais para todos os níveis de agregação. A equação para a abordagem *top-down* será, então:

$$
\mathbfit{\tilde{y}}_{T+h | T} = \mathbfit{SG\hat{y}}_{T+h | T}
$$ {#eq-topdown_3}

Na notação matricial para a estrutura da @fig-h, temos:

$$
\begin{bmatrix}
    \tilde{y}_{t} \\
    \tilde{y}_{A, t} \\
    \tilde{y}_{B, t} \\
    \tilde{y}_{C, t} \\
    \tilde{y}_{AA, t} \\
    \tilde{y}_{AB, t} \\
    \tilde{y}_{BA, t} \\
    \tilde{y}_{BB, t} \\
    \tilde{y}_{CA, t} \\
    \tilde{y}_{CB, t}
\end{bmatrix}_{10 \times 1}
=
\mathbfit{S}_{10 \times 6}
\begin{bmatrix}
    p_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}_{6 \times 10}
\begin{bmatrix}
    \hat{y}_{T+h|T} \\
    \hat{y}_{A, T+h|T} \\
    \hat{y}_{B, T+h|T} \\
    \hat{y}_{C, T+h|T} \\
    \hat{y}_{AA, T+h|T} \\
    \hat{y}_{AB, T+h|T} \\
    \hat{y}_{BA, T+h|T} \\
    \hat{y}_{BB, T+h|T} \\
    \hat{y}_{CA, T+h|T} \\
    \hat{y}_{CB, T+h|T}
\end{bmatrix}_{10 \times 1}
$$ {#eq-matriz_topdown1}

Multiplicando as matrizes $\mathbfit{G}$ e $\mathbfit{\hat{y}}$ obtemos o vetor $\mathbfit{b}_t$^[Ver @eq-vetor_b.], que contém as previsões para os níveis mais desagregados, agora definidas como proporção do total.
$$
\begin{bmatrix}
    \tilde{y}_{t} \\
    \tilde{y}_{A, t} \\
    \tilde{y}_{B, t} \\
    \tilde{y}_{C, t} \\
    \tilde{y}_{AA, t} \\
    \tilde{y}_{AB, t} \\
    \tilde{y}_{BA, t} \\
    \tilde{y}_{BB, t} \\
    \tilde{y}_{CA, t} \\
    \tilde{y}_{CB, t}
\end{bmatrix}_{10 \times 1}
=
\mathbfit{S}_{10 \times 6}
\begin{bmatrix}
    p_1\hat{y}_{T+h|T} \\
    p_2\hat{y}_{T+h|T} \\
    p_3\hat{y}_{T+h|T} \\
    p_4\hat{y}_{T+h|T} \\
    p_5\hat{y}_{T+h|T} \\
    p_6\hat{y}_{T+h|T}
\end{bmatrix}_{6 \times 1}
$$ {#eq-matriz_topdown2}

Note que, por se tratar de um método *top-down*, é necessário apenas o primeiro elemento do vetor de previsões base, ou seja, a previsão do nível mais agregado (@eq-matriz_topdown2). Sendo essa exatamente uma das vantagens do método *top-down*, na prática, podemos anular os demais elementos de $\mathbfit{\hat{y}}$.

$$
\begin{bmatrix}
    \tilde{y}_{t} \\
    \tilde{y}_{A, t} \\
    \tilde{y}_{B, t} \\
    \tilde{y}_{C, t} \\
    \tilde{y}_{AA, t} \\
    \tilde{y}_{AB, t} \\
    \tilde{y}_{BA, t} \\
    \tilde{y}_{BB, t} \\
    \tilde{y}_{CA, t} \\
    \tilde{y}_{CB, t}
\end{bmatrix}_{10 \times 1}
=
\mathbfit{S}_{10 \times 6}
\begin{bmatrix}
    p_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    p_6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}_{6 \times 10}
\begin{bmatrix}
    \hat{y}_{T+h|T} \\
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    0 \\
    0
\end{bmatrix}_{10 \times 1}
$$ {#eq-matriz_topdown3}

Substituindo a matriz $\mathbfit{S}$, temos as equações que definem as previsões reconciliadas.

$$
\begin{bmatrix}
    \tilde{y}_{t} \\
    \tilde{y}_{A, t} \\
    \tilde{y}_{B, t} \\
    \tilde{y}_{C, t} \\
    \tilde{y}_{AA, t} \\
    \tilde{y}_{AB, t} \\
    \tilde{y}_{BA, t} \\
    \tilde{y}_{BB, t} \\
    \tilde{y}_{CA, t} \\
    \tilde{y}_{CB, t}
\end{bmatrix}_{10 \times 1}
=
\begin{bmatrix}
    1 & 1 & 1 & 1 & 1 & 1 \\
    1 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 1 \\
    1 & 0 & 0 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}_{10 \times 6}
\begin{bmatrix}
    p_1\hat{y}_{T+h|T} \\
    p_2\hat{y}_{T+h|T} \\
    p_3\hat{y}_{T+h|T} \\
    p_4\hat{y}_{T+h|T} \\
    p_5\hat{y}_{T+h|T} \\
    p_6\hat{y}_{T+h|T}
\end{bmatrix}_{6 \times 1}
$$ {#eq-matriz_topdown3}

Já a abordagem *bottom-up* parte do raciocínio inverso e define as previsões de cada elemento da estrutura a partir das previsões dos elementos mais desagregados. Para tanto, basta modificar a matriz $\mathbfit{G}$.

$$
\mathbfit{G}
=
\begin{bmatrix}
    0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}_{6 \times 10}
$$ {#eq-matriz_gbu}

Portanto, $\mathbfit{G}$ define a abordagem --- se *top-down*, *bottom-up* ou outro qualquer ---, e $\mathbfit{S}$ define a maneira da qual as previsões são somadas para formar as equações de previsão para cada elemento da estrutura. Portanto, chamamos $\mathbfit{G}$ de matriz de reconciliação.

$$
\begin{bmatrix}
    \tilde{y}_{t} \\
    \tilde{y}_{A, t} \\
    \tilde{y}_{B, t} \\
    \tilde{y}_{C, t} \\
    \tilde{y}_{AA, t} \\
    \tilde{y}_{AB, t} \\
    \tilde{y}_{BA, t} \\
    \tilde{y}_{BB, t} \\
    \tilde{y}_{CA, t} \\
    \tilde{y}_{CB, t}
\end{bmatrix}_{10 \times 1}
=
\begin{bmatrix}
    1 & 1 & 1 & 1 & 1 & 1 \\
    1 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 1 \\
    1 & 0 & 0 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}_{10 \times 6}
\begin{bmatrix}
    \hat{y}_{AA, T+h|T} \\
    \hat{y}_{AB, T+h|T} \\
    \hat{y}_{BA, T+h|T} \\
    \hat{y}_{BB, T+h|T} \\
    \hat{y}_{CA, T+h|T} \\
    \hat{y}_{CB, T+h|T}
\end{bmatrix}_{6 \times 1}
$$ {#eq-matriz_bottomup}

Quando $m$ — a quantidade de elementos do nível mais desagregado — é muito grande, tornando muito custoso obter $\mathbfit{\hat{y}_t}$, e não se deseja uma abordagem estritamente *top-down*, pode-se combinar as duas formas:

$$
\mathbfit{G}
=
\begin{bmatrix}
    0 & p_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    0 & p_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    0 & 0 & p_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    0 & 0 & p_4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
    0 & 0 & 0 & p_5 & 0 & 0 & 0 & 0 & 0 & 0 \\
    0 & 0 & 0 & p_6 & 0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}_{m \times n}
$$ {#eq-matriz_gmo}

Esse método é chamado de *middle-out*. Nele, o vetor $\mathbfit{b}_t$ é reescrito como proporções de algum nível intermédiário arbitrariamente escolhido, ao invés de proporções do total. Isso permite uma abordagem mais econômica, em termos de custo computacional e de tempo, ao mesmo tempo em que mantém em algum grau as características individuais das hierarquias.

$$
\begin{bmatrix}
    \tilde{y}_{t} \\
    \tilde{y}_{A, t} \\
    \tilde{y}_{B, t} \\
    \tilde{y}_{C, t} \\
    \tilde{y}_{AA, t} \\
    \tilde{y}_{AB, t} \\
    \tilde{y}_{BA, t} \\
    \tilde{y}_{BB, t} \\
    \tilde{y}_{CA, t} \\
    \tilde{y}_{CB, t}
\end{bmatrix}_{n \times 1}
=
\begin{bmatrix}
    1 & 1 & 1 & 1 & 1 & 1 \\
    1 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 1 \\
    1 & 0 & 0 & 0 & 0 & 0 \\
    0 & 1 & 0 & 0 & 0 & 0 \\
    0 & 0 & 1 & 0 & 0 & 0 \\
    0 & 0 & 0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 0 & 1 & 0 \\
    0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}_{n \times m}
\begin{bmatrix}
    p_1\hat{y}_{A,T+h|T} \\
    p_2\hat{y}_{A,T+h|T} \\
    p_3\hat{y}_{A,T+h|T} \\
    p_4\hat{y}_{B,T+h|T} \\
    p_5\hat{y}_{B,T+h|T} \\
    p_6\hat{y}_{B,T+h|T}
\end{bmatrix}_{m \times 1}
$$ {#eq-matriz_mo}

## Coerência e reconciliação

Seja somando as previsões do nível mais desagregado para formar os níveis superiores da hierarquia (*bottom-up*) ou distribuindo proporcionalmente as previsões do nível mais agregado (*top-down*), o vetor $\mathbfit{\tilde{y}}_t$ representa as previsões *coerentes*. Isso significa que as previsões são totalizadas corretamente --- as previsões de cada elemento agregado corresponde ao somatório das previsões dos níveis inferiores da hierarquia. Isso é garantido pela multiplicação das matrizes $\mathbfit{SG}$.

Não fosse essa pré-multiplicação, nada garantiria a coerência das previsões. Tomando a estrutura da @fig-h como exemplo, seria um acaso improvável que as previsões do agregado para o estado do Espírito Santo fossem exatamente a soma das previsões individuais de seus municípios. Isso porque não há qualquer razão para que cada série siga o mesmo processo estocástico.

Os métodos para gerar previsões coerentes ($\mathbfit{\tilde{y}_t}$) a partir de previsões base ($\mathbfit{\hat{y}_t}$) são chamados de métodos de *reconciliação*. Os métodos de reconciliação tradicionais apresentados, *top-down* e *bottom-up*, utilizam informação limitada. No método *top-down*, utiliza-se apenas informações do nível mais agregado --- por isso, apenas a primeira coluna em (@eq-matriz_g) é diferente de zero. Já na abordagem *bottom-up*, utiliza-se apenas as informações dos níveis mais desagregados, o que resulta na submatriz identidade $m \times m$ (@eq-matriz_gbu), enquanto as colunas que representam os níveis mais agregados são nulas.

Alternativamente, podemos pensar numa matriz $\mathbfit{G}$ qualquer que utilize toda a informação disponível e tenha algumas propriedades que garantam que as previsões coerentes tenham o menor erro o possível. Esse é o problema de pesquisa trabalhado na *reconciliação ótima*.

Os métodos analíticos de reconciliação ótima são aqueles que estimam a matriz de reconciliação $\mathbfit{G}$ através de regressão linear. Isso resulta na redefinição das previsões do nível mais desagregado como uma combinação linear^[Por essa razão, esses métodos são também chamados de métodos de combinação.] das previsões de todos os elementos de todos os níveis da hierarquia, utilizando, assim, toda a informação disponível.

A primeira abordagem prática para reconciliação ótima, via MQO, é formalizada apenas em @hyndman_optimal_2011. @hyndman_fast_2016 tentam aperfeiçoar o método usando as variâncias das previsões individuais estimadas (dentro da amostra) como estimativa para a matriz de variância-covariância dos erros de reconciliação, de forma a utilizá-las como pesos e realizar a reconciliação ótima por mínimos quadrados ponderados (MQP).

@wickramasuriya_optimal_2019 argumentam que o que de fato interessa é que as previsões reconciliadas tenham o menor erro. Então, corrigem a abordagem de reconciliação ótima para o objetivo de minimização dos erros das previsões reconciliadas $\mathbfit{\tilde{y}_{t+h}}$, ao invés dos erros das previsões individuais $\mathbfit{\hat{y}_{t+h}}$. Dado que isso implica na minimização da variância de $\mathbfit{\tilde{e}_{t+h}}$, ou seja, na minimização do somatório da diagonal, o traço, da matriz de variância-covariância de $\mathbfit{\tilde{e}_{t+h}}$, eles chamaram esse método de Traço Mínimo (MinT, na sigla em inglês). 

Este método é considerado o estado-da-arte para a reconciliação ótima analítica. Nele, o objetivo é minimizar o erro das previsões reconciliadas:

$$
\mathbfit{\tilde{e}}_{t+h|t} = \mathbfit{y}_{t+h} - \mathbfit{\tilde{y}}_{t+h|t}
$$ {#eq-mint1}

Essa equação pode ser reescrita como $\mathbfit{\tilde{e}}_t = \mathbfit{SG\hat{e}}_t$, que tem variância dada por^[Ver demonstração \ref{proposicao3}.]

$$
\text{Var}[\mathbfit{\tilde{e}}] = \mathbfit{SG\hat{W}}_{t+h|t}\mathbfit{G'S'}
$$ {#eq-mint2}

\noindent em que $\mathbfit{\hat{W}}_{t+h|t}$ é a matriz de variância-covariância dos erros de previsão base.

A abordagem consiste então em se obter um valor de $\mathbfit{G}$ que minimize o traço de $\text{Var}[\mathbfit{y}_{t+h} - \mathbfit{\tilde{y}}_{t+h|t}]$. Isso resultaria no melhor (variância mínima) estimador linear não viesado^[A ausência de viés é garantida pela \nameref{proposicao1}.] [@wickramasuriya_optimal_2019]. Sua solução é dada por

$$
\mathbfit{G} = (\mathbfit{S'\hat{W}}^\dagger_{t+h|t}\mathbfit{S})^{-1}\mathbfit{S'\hat{W}}^\dagger_{t+h|t}
$$ {#eq-mint3}

\noindent em que $\mathbfit{\hat{W}}^\dagger_{t+h|t}$ é a inversa generalizada de Moore-Penrose para $\mathbfit{\hat{W}}_{t+h|t}$^[A necessidade da inversa generalizada aqui é trivial, uma vez que a inversa regular, do tipo $\mathbfit{A}^{-1}$, requer matriz quadrada, o que não acontece no caso de séries temporais hierárquicas, uma vez que, necessariamente, tem-se $n>m$ ($n=m+\text{número de nós de agregação}$). Além disso, $\mathbfit{\hat{W}}$ é posto incompleto (ver demonstração \ref{proposicao4}).]. Essa formulação corresponde a um problema de regressão por mínimos quadradados generalizados, considerando $\mathbfit{S}$ como a matriz de preditores e $\mathbfit{G}$ os coeficientes a serem estimados. Consequentemente, as previsões ótimas reconciliadas são dadas por

$$
\mathbfit{\tilde{y}}_{t+h|t} = \mathbfit{S}(\mathbfit{S'\hat{W}}^\dagger_{t+h|t}\mathbfit{S})^{-1}\mathbfit{S'\hat{W}}^\dagger_{t+h|t}\mathbfit{\hat{y}}_{t+h|t}
$$ {#eq-mint4}

\noindent que são ao menos tão boas quanto as previsões individuais, de acordo com o demonstrado pelos autores.

@spiliotis_hierarchical_2021 propõem a utilização de *machine learning* para a reconciliação ótima de séries temporais, especificamente os métodos de RF e *gradient boosting*. A ideia subjacente é treinar modelos individuais de ML para cada série do nível mais desagregado a partir de previsões para dentro da amostra obtidas via *rolling forecasting*, com o objetivo de encontrar a combinação contemporânea que minimiza os erros, para então usar esses modelos para combinar as previsões base. Os autores descrevem como vantagens desse método em relação aos anteriores a descrição de relacionamentos não lineares, performance preditiva e a desnecessidade da utilização de todos os elementos da hierarquia na combinação ótima. Para o conjunto de dados utilizados, os autores afirmam que o uso de ML, especialmente o XGB, permite alcançar maior performance em relação às abordagens de nível único e o MinT.

@panagiotelis_forecast_2021 reinterpretam a literatura de coerência e reconciliação de previsões pontuais a partir de uma abordagem geométrica, trazendo provas alternativas para conclusões anteriores ao mesmo tempo em que fornece novos teoremas. Além disso, os autores estendem essa interpretação geométrica para o contexto probabilístico, fornecendo métodos paramétricos e não paramétricos (via *bootstrapping*) para reconciliação de previsões probabilísticas.

@shiratori_prediction_2020 propõem a utilização de um termo de regularização para penalizar o erro dos níveis superiores na previsão das séries do nível mais desagregado, de forma a obter as previsões base já com a informação dos níveis superiores, substituindo a abordagem de combinação ótima. Assim como nas demais abordagens, a coerência é obtida através de agregação BU.

Uma revisão sistemática descrevendo os avanços e diferentes abordagens para reconciliação de previsões hierárquicas pode ser encontrada em @athanasopoulos_forecast_2023. Isto posto, a presente dissertação foca em estender a abordagem de combinação ótima por ML, realizada em trabalhos como de @spiliotis_hierarchical_2021, para a reconciliação de séries temporais hierárquicas e agrupadas, avaliando a efetividade de 7 implementações de métodos de ML (RF, XGB, *Lasso*, *ridge*, *elastic net*, SVM e LGBM), propondo variações metodológicas para a obtenção do conjunto de treinamento dos modelos de ML e obtendo ganhos significativos de performance em três estudos de caso.

# ESTENDENDO O *FRAMEWORK* DE COMBINAÇÃO ÓTIMA VIA *MACHINE LEARNING*

## Introdução

A ideia da reconciliação de previsões por combinação consiste em encontrar combinações, lineares ou não, que minimizem o erro de combinação contemporânea dentro da amostra e tenham o menor erro de generalização para o conjunto de previsões base, incoerentes, para fora da amostra. O objetivo é que, dada a performance das previsões base individuais $\hat{y}_t$, as previsões reconciliadas $\tilde{y}_t$, obtidas por combinação e agregadas de maneira BU, se tornem mais precisas por conta da incorporação de informações dos diferentes níveis hierárquicos e demais nós da estrutura. 

Embora os métodos de ML também sejam utilizados no contexto de previsão de séries temporais, é importante ressaltar que este não é caso deste trabalho. As previsões de séries temporais são realizadas anteriormente, na obtenção das previsões base $\hat{y}_t$. Nada impede que essas previsões sejam obtidas por métodos de ML, inclusive utilizando técnicas de previsão hierárquica na predição de séries individuais [@athanasopoulos_forecasting_2017], entretanto, a aplicação dos métodos de ML aqui ocorrem na reconciliação, ou seja, na combinação contemporânea, *cross section*, das previsões base que, por sua vez, podem ter sido obtidas através de qualquer método.

Nesse sentido, este trabalho busca estender o *framework* de combinação ótima utilizando métodos de ML para a reconciliação de séries temporais hierárquicas e agrupadas. Além de avaliar 7 implementações de métodos de ML (RF, XGB, *Lasso*, *ridge*, *elastic net*, SVM e LGBM)^[Uma breve introdução sobre cada um desses métodos pode ser encontrado no Apêndice \ref{apendice_metodos_ml}.], propõe variações metodológicas para a obtenção do conjunto de treinamento dos modelos de ML.

## Metodologia de combinação ótima via *machine learning* {#sec-metodologia_ml}

@spiliotis_hierarchical_2021 propõem estimar essa combinação em uma abordagem que constitui-se nas seguintes etapas:

@. Previsões base:
Para cada série $y_t$, incluindo os nós de agregação, utiliza-se todo o conjunto de treinamento, de tamanho $T$, para produzir previsões $h$ passos à frente, $\hat{y}_{t+h|T}$, utilizando qualquer método de previsão. Essas são as previsões base a serem reconciliadas.

![Previsões base a serem reconciliadas. Em azul, o conjunto de treino e, em verde, as previsões $h$ passos à frente para fora da amostra.](img/modelagem_1.png){#fig-modelagem-1 fig-scap="Previsões base a serem reconciliadas" width=70%}

@. *Rolling forecast*:
A amostra treino é separada em $p$ subamostras para todas as séries em todos os níveis de agregação, com a primeira separação abrangindo $y_1$ até $y_Q$, a segunda $y_1$ até $y_{Q+1}$ e assim sucessivamente até a última separação abrangendo $y_1$ até $y_{Q+p-1}$. Para cada subamostra é treinado um modelo que, mais uma vez, pode ser estimado utilizando qualquer método. Então, é realizada previsão um passo a frente. Essas previsões são então concatenadas para formar o conjunto de treinamento para os modelos de ML.

![Esquema de modelagem *rolling forecast*](img/modelagem_2.png){#fig-modelagem-2 width=70%}

@. Treino dos modelos de ML:
Para cada série do nível mais desagregado, $y_m$, é treinado um modelo de ML com $n+1$ variáveis, compostas pelas $n$ séries previstas — que incluem todos os níveis de agregação —, mais a própria $y_m$ como variável explicada (@tbl-modelagem). Cada uma das $n$ séries contam com $p$ previsões obtidas na 2ª etapa, representadas na @fig-modelagem-2 pela cor verde. Isso resulta em um modelo de reconciliação ótima para cada elemento do menor nível da hierarquia, combinando informações disponíveis de todos os níveis hierárquicos.

```{r}
#| label: tbl-modelagem
#| tbl-cap: "Conjunto de treinamento para predição dos modelos de ML"
  tibble::tribble(
    ~`Var. Explicada`, ~`Variável 1`, ~`Variável 2`, ~..., ~`Variável n`,
    "$y_{1,Q+1}$", "$\\hat{y}_{1,Q+1}$", "$\\hat{y}_{2,Q+1}$", "...", "$\\hat{y}_{n,Q+1}$",
    "$y_{1,Q+2}$", "$\\hat{y}_{1,Q+2}$", "$\\hat{y}_{2,Q+2}$", "...", "$\\hat{y}_{n,Q+2}$",
    "...", "...", "...", "...", "...",
    "$y_{m,Q+p}$", "$\\hat{y}_{1,Q+p}$", "$\\hat{y}_{2,Q+p}$", "...", "$\\hat{y}_{n,Q+p}$"
  ) |>
  kbl(booktabs = TRUE, escape = FALSE) |>
  kable_styling(latex_options = c("striped"))
```

@. Reconciliação ótima:
Com os modelos treinados na 3ª etapa, passa-se as previsões base obtidas na 1ª etapa como regressores para se obter as previsões reconciliadas $h$ passos à frente das séries do nível mais desagregado $\tilde{y}_m$.

@. Agregação: Assim como nos métodos analíticos de combinação ótima, a obtenção das previsões reconciliadas para os demais níveis de hierárquicos $\tilde{y}_n$ se dá através da agregação BU, mas ao invés de se somar as previsões base $\hat{y}_m$, somam-se as previsões reconciliadas $\tilde{y}_m$.

Dessa forma, essa metodologia é semelhante à aplicada na reconciliação ótima analítica, se afastando principalmente em três pontos: (i) a utilização de algoritmos de ML ao invés de MQG, (ii) a não atribuição de peso de forma obrigatória para todos os nós da hierarquia e (iii) no ajuste de um modelo individual para cada série do nível mais desagregado, permitindo maior especialização e sendo capaz de se adaptar melhor aos diferentes padrões de cada série [@spiliotis_hierarchical_2021]. A @fig-modelagem-rolling ilustra a metodologia de combinação ótima via ML.

![Metodologia de combinação ótima via ML. Estratégia *rolling forecast* (adaptado de @spiliotis_hierarchical_2021)](img/modelagem_rolling.png){#fig-modelagem-rolling fig-scap="Metodologia de combinação ótima via ML. Estratégia *rolling forecast*"}

Um ponto negativo na estratégia de *rolling forecasting* (2ª etapa) é o custo em tempo de processamento. Em um cenário de hierarquias largas e agrupadas, o ganho de performance com a estimação de um modelo individual para cada $p$ subamostras pode não compensar o *trade-off* em relação ao custo.

Nesse sentido, uma alternativa é o processo de reestimação, que consiste na reestimação de um modelo para um novo conjunto de dados. Nessa estratégia, que chamarei de *reduced fitted values*, os $n$ modelos são treinados com o conjunto de dados $\Omega_1=\{y_1, y_2, ..., y_Q\}$ e então reajustados para o restante do conjunto de dados $\Omega_2=\{y_{Q+1}, y_{Q+2}, ..., y_{Q+p}\}$, se esquivando do processo de *rolling origin*. De posse dos modelos reajustados, utiliza-se então os valores ajustados (*fitted values*) como conjunto de treinamento para os modelos de ML. Uma restrição dessa abordagem é que, fixados os hiperparâmetros de $\hat{y}_{t|\Omega_1}$, não necessariamente todos os coeficientes de $\hat{y}_{t|\Omega_2}$ alcançarão convergência no reajuste de seus coeficientes para o novo conjunto de dados.

A metodologia de combinação ótima na estratégia *reduced fitted values* está ilustrada na @fig-modelagem-refit. Nela, o tamanho da amostra disponível para o treino dos modelos de ML é o mesmo que o da estratégia por *rolling forecasting*, mas o custo computacional é reduzido, uma vez que é treinado apenas um modelo para cada série $y_t$. Além disso, como alterações são comuns em processos estocásticos de séries longas de natureza econômica ou social, a estratégia *reduced fitted values* pode resultar em previsões mais precisas, uma vez que os modelos são reestimados apenas com os dados de $\Omega_2$, resultando em coeficientes adaptados às características mais recentes das séries.

![Metodologia de combinação ótima via ML. Estratégia *reduced fitted values*](img/modelagem_refit.png){#fig-modelagem-refit}

As duas estratégias possuem uma dificuldade em comum: a necessidade de realização de previsões para dentro do conjunto de treinamento. Isso pode ser um problema para séries temporais com poucas observações ou incompletas. Por exemplo, suponha o problema de previsão hierárquica para uma franquia de varejo. O fechamento e abertura de novas filiais é bastante frequente, o que eventualmente irá resultar em filiais abertas após o período escolhido para $Q$. Como a combinação ótima trata da combinação contemporânea, isso significaria lidar com dados faltantes no conjunto de treinamento dos modelos de ML. As alternativas metodológicas para contornar esse problema passam por restringir a escolha dos métodos de ML para apenas aqueles capazes de trabalhar com dados faltantes nativamente ou utilizar uma estratégia de imputação na etapa de pré-processamento.

Para as séries de processo estocástico bem definido, mesmo em longos períodos, uma 3º estratégia para aquisição do conjunto de treinamento para os modelos de ML é usar os valores ajustados dos modelos das previsões base $\hat{y}$ (passo 1) como input para os modelos de ML. Chamarei essa estratégia de *fitted values* (@fig-modelagem-fitted).

![Metodologia de combinação ótima via ML. Estratégia *fitted values*](img/modelagem_fitted.png){#fig-modelagem-fitted}

Essa estratégia equivale a usar previsões um passo à frente obtidas a partir do conjunto de dados completo $\Omega = \Omega_1 + \Omega_2$. Em relação às duas estratégias anteriores, essa abordagem aumenta o tamanho da amostra treino em $Q$ observações para cada série, o que pode melhorar a qualidade modelos de ML. Além disso, ela é ainda mais simples de implementar, uma vez que não é necessário gerar novos modelos, utilizando aqueles já estimados na 1ª etapa. 

## Notas sobre otimização de hiperparâmetros na reconciliação ótima

Considere uma função de ajuste $f$, um conjunto de pontos $\Omega = {\omega_1, ..., \omega_n}$ com $\omega_i = (\mathbf{x}_i \; \mathbf{y}_i)'$, variáveis de decisão ou parâmetros $\mathbf{x}_i \in \mathbb{R}^m$ e imagem $\mathbf{y}_i = f(\mathbf{x}_i) \in \mathbb{R}$. Diferentemente da abordagem do modelo clássico de regressão linear em que, cumpridas certas hipóteses, há um modelo teórico de coeficientes estimados por mínimos quadrados ordinários (MQO) que é garantido pelo teorema de Gauss-Markov ser o melhor estimador linear não viesado (BLUE), em ML o objetivo é encontrar, de forma iterativa, um modelo que melhor aproxima a função $f$ usando a informação contida em $\Omega$, ou seja, queremos ajustar uma função de regressão $\hat{f}_\Omega$ aos nossos dados $\Omega$ de forma que $\hat{\mathbf{y}} = \hat{f}_\Omega(\mathbf{x}, \varepsilon)$ tenha o menor erro de aproximação $\varepsilon$ [@bischl_resampling_2012].

Para verificar o quão bem o modelo $\hat{f}_D$ se aproxima da função real $f$, é necessário uma função de perda $L(\mathbf{y}, \hat{f}(\mathbf{x}))$ que, no caso de regressão, será a perda quadrática $(\mathbf{y} - \hat{f}(\mathbf{x}))^2$ ou a perda absoluta $|\mathbf{\mathbf{x}} - \hat{f}(\mathbf{x})|$. Esses valores são agregados pela média para formar as funções de custo erro médio quadrático (MSE) e erro médio absoluto (MAE).

A fim de evitar o problema de sobreajuste, diferentes métodos para mensuração do erro de generalização (i.e. performance fora da amostra) foram desenvolvidos, sendo a validação cruzada [@stone_cross-validation_1974] um dos mais populares. Este método consiste em gerar repetidamente $i$ subconjuntos de treino $\Omega_i^a$ e teste $\Omega_i^b$, com $\Omega = \Omega_i^a + \Omega_i^b$; ajustar um meta-modelo com cada conjunto de treino; e avaliar sua qualidade no conjunto de teste correspondente. Dessa forma, a estimativa do erro de generalização se dá por

$$
\widehat{GE}_{me} = \frac{1}{k}\sum_{i=1}^{k}\widehat{GE}(\hat{f}_{\Omega_i^a}, \Omega_i^b)
$$ {#eq-ge-cv}

\noindent sendo $\widehat{GE}_{me}$ o erro de generalização médio. Dividindo a amostra em $k$ subconjuntos, utilizando $k-1$ para ajustar um meta-modelo e validando a performance no subconjunto restante — e repetindo esse procedimento para todas as possibilidades de subconjuntos —, temos a validação cruzada $k$-*fold*. Esta é uma técnica de reamostragem que permite estimar o erro de generalização de um modelo de forma mais robusta e confiável que o *hold-out* simples (quando $k=1$). Isso porque, ao contrário do *hold-out*, a validação cruzada permite que todos os dados sejam usados tanto para treino quanto para teste, o que reduz a variância da estimativa do erro de generalização.

Para dados *cross-section*, essa definição de validação cruzada $k$-*fold* é geralmente suficiente. Entretanto, no caso de séries temporais, o analista deve tomar alguns cuidados na escolha da abordagem de validação cruzada. Isso porque, ao contrário de dados *cross-section*, os dados de séries temporais são dependentes no tempo. Isso traz dois problemas aparentes: primeiramente, ao dividir a amostra em $k$ subconjuntos aleatórios, o meta-modelo será treinado em períodos descontínuos e com dados futuros aos dados de teste em, ao menos, $k-1$ subconjuntos. Em segundo lugar, realizando o subdivisão treino-teste em $y_t$, com $y_{t-1}$ na amostra treino, significa que os subconjuntos de treino e teste são dependentes (@fig-cv-aleatoria).

![Validação $k$-fold aleatória](img/cv-aleatorio.png){#fig-cv-aleatoria width=70%}

Para contornar esse problema, outros métodos de validação cruzada foram desenvolvidos pensando em séries temporais. Tomando o conjunto de validação com dados exclusivamente posteriores aos dados de treino, temos o método conhecido como validação cruzada com origem móvel (@fig-cv-blocked).

![Validação $k$-fold com origem móvel](img/cv-blocked.png){#fig-cv-blocked width=70%}

Já quando excluímos as defasagens dependentes do subconjunto de treino (em relação ao conjunto de validação), temos a validação cruzada não dependente (@fig-cv-n-dep). Para ilustrar esse caso, tome um processo AR(3). Excluir as defasagens dependentes significa que, se o conjunto de validação começa em $y_t$, então o conjunto de treino pode conter observações apenas até $y_{t-4}$. O problema evidente dessa abordagem é que, dependendo do tamanho da estrutura de autocorrelação da série, muitas defasagens são excluídas de cada *k-fold*, podendo inviabilizar o processo de validação.

![Validação $k$-fold não-dependente](img/cv-n-dep.png){#fig-cv-n-dep width=70%}

Cumpridas algumas restrições, esses problemas são apenas aparentes. @bergmeir_use_2012 comparam uma série de métodos de validação em séries temporais e concluem que, considerando séries estacionárias, os problemas teóricos relacionados à dependência não produzem impactos empíricos. @bergmeir_note_2018 vão além e abordam esse problema teórica e empiricamente, concluindo que não apenas é possível o uso de validação cruzada $k$-*fold* quando a série não apresenta autocorrelação serial nos resíduos, mas também é uma melhor escolha do que a validação fora da amostra^[*Out-of-sample evaluation* é o método padrão na literatura de séries temporais e consiste na separação da porção final da série — geralmente entre 20% e 30% — para validação.].

Entretanto, no mundo real, trabalhando com uma longa coleção de séries temporais hierárquicas agrupadas, raramente será o caso em que a totalidade das séries são estacionárias. Nesse sentido, é razoável supor que, mesmo que o método de reconciliação via combinação ótima trate da combinação contemporânea das previsões, se as médias das séries estão se alterando ao longo do tempo, a combinação ótima entre elas também possa se alterar. Por isso, um processo de validação cruzada que não suponha estacionariedade pode ser mais adequado.

Por essa razão, proponho que o processo de reamostragem para otimização de hiperparâmetros conte com, além da validação cruzada $k$-fold, a validação cruzada com origem móvel em seu *benchmark*. Para tanto, é necessário definir uma janela mínima para $\Omega^a$ e um horizonte para o conjunto de teste $\Omega^b$, que devem ser escolhidos de acordo com as características do conjunto de dados trabalhado.

## Conclusão

Nesta seção, foram apresentadas a abordagem para a reconciliação de previsões de séries temporais hierárquicas e agrupadas por combinação ótima utilizando métodos de ML. Foram propostas três estratégias para a obtenção do conjunto de treinamento dos modelos de ML: *rolling forecast*, *reduced fitted values* e *fitted values*. Cada uma dessas estratégias tem suas vantagens e desvantagens, e a escolha entre elas depende das características específicas das séries temporais que estão sendo previstas.

Além disso, discutiu-se a importância da validação cruzada na otimização de hiperparâmetros para modelos de ML no cenário de combinação ótima. Em particular, destacou-se a necessidade de considerar a dependência temporal nos dados ao escolher o método de validação cruzada. Nesse sentido, propus o uso da validação cruzada com origem móvel como uma alternativa para calibração dos modelos de combinação ótima.

# PREVISÕES DE SALDOS DE EMPRÉSTIMOS E FINANCIAMENTOS DO BANCO DO ESTADO DO ESPÍRITO SANTO {#sec-estban}

## Introdução

Embora no séc. XX ainda houvesse espaço para uma gestão guiada apenas por instinto [@wallander_budgeting_1999], atualmente é impensável um banco não realizar previsões de seus resultados e comunicar suas expectativas ao mercado. Nesse documento, ou *guidance*, a projeção da carteira de crédito — o total de empréstimos e financiamentos, dentre outros itens — é frequentemente a primeira informação fornecida. Juntamente com as projeções de depósitos, provisões para créditos de liquidação duvidosa, eficiência operacional, entre outros indicadores-chave, essas projeções determinam a temperatura das expectativas da instituição, e isso é essencial para os acionistas e investidores. Essas projeções precisam ser tão precisas quanto possível para que se possa calcular o risco de transacionar com a instituição financeira.

Ainda que não existam penalidades específicas para instituições financeiras que erram (por uma boa margem) em suas projeções, elas podem sofrer consequências negativas em outros aspectos, como na avaliação de seus desempenhos por parte dos investidores e clientes. Estes podem considerar as projeções equivocadas como um sinal de falta de competência ou confiança na instituição financeira, o que pode afetar negativamente a reputação e a imagem da instituição.

Além disso, nos casos em que algum grupo se sentir lesado, os bancos podem enfrentar ações judiciais se suas projeções forem consideradas enganosas ou fraudulentas. Por exemplo, se uma instituição financeira fizer projeções excessivamente otimistas para incentivar os investidores a comprar seus títulos e, posteriormente, as projeções se mostrarem incorretas, ela pode ser acusada de fraude^[Art. 3º: Divulgar informação falsa ou prejudicialmente incompleta sobre instituição financeira. Pena: Reclusão, de 2 (dois) a 6 (seis) anos, e multa. Art. 4º: Gerir fraudulentamente instituição financeira. Pena: Reclusão, de 3 (três) a 12 (doze) anos, e multa [@brasil_lei_1986].] ou, ao menos, gestão temerária^[Art. 4º, parágrafo único: Se a gestão é temerária: Pena: Reclusão, de 2 (dois) a 8 (oito) anos, e multa [@brasil_lei_1986].] — ambos caracterizados como crime contra o Sistema Financeiro Nacional (SFN).

Por isso, é importante que as instituições financeiras sejam transparentes e precisas em suas projeções, fornecendo informações confiáveis e atualizadas para seus clientes e investidores. No entanto, há também motivações estratégicas para essa atividade. @beccalli_earnings_2015 mostraram que, em uma amostra de 55 bancos europeus, a utilização de *guidance* está associada a um aumento de 15% na probabilidade do banco atingir ou superar as expectativas de mercado. Isso, por sua vez, está associado a um incremento de até 5% no retorno por ação em relação aos bancos que não alcançaram ou superaram as expectativas.

A prática usual em *budgeting*^[O orçamento é um documento no qual é definido o planejamento financeiro de um empresa, geralmente para o ano seguinte, estabelecendo metas e objetivos. Nele, são projetadas as expectativas da empresa e é base de comparação para saber como os resultados estão se desviando da performance esperada.], principalmente para empresas com muitas filiais, é a *top-down*. No caso dos bancos de varejo, com muitas agências espalhadas pelo território, esse método pode ser muito útil. Na prática, em um país continental como o Brasil, a quantidade de filiais para as quais se deve realizar previsões pode tornar a tarefa desafiadora mesmo para instituições regionais.

Esse é o caso do Banestes. Com 134 agências distribuídas pelos 78 municípios capixabas, realizar o *budgeting* para R$ 5,5 bi de faturamento^[Conforme demonstrativos publicados referentes ao exercício de 2022 [@banco_do_estado_do_espirito_santo_demonstracoes_2022].] não é uma tarefa trivial. Além de uma estrutura hierárquica de alta dimensionalidade por conta da quantidade de agências, se tratando de um banco múltiplo^[Para ser classificado como banco múltiplo, a instituição financeira deve operar com, no mínimo, duas carteiras dentre: comercial; investimento ou desenvolvimento; crédito imobiliário; de crédito, financiamento e investimento, e; arrendamento mercantil [@conselho_monetario_nacional_resolucao_1994].] que opera com diversas carteiras, as $n$ modalidades de crédito^[Crédito consignado, rural, imobiliário, pessoal, capital de giro, desconto de títulos etc.] expandem a estrutura para um total de $n \times 134$ séries temporais a serem estimadas.

Dada tal complexidade, a abordagem *top-down* se coloca como uma opção viável em termos de tempo de processamento e análise. No entanto, conforme descemos na hierarquia, menos precisa ela se torna e, além disso, as características individuais das séries temporais do menor nível hierárquicos são ignoradas. No sentido inverso, caso se escolha estimar modelos para cada uma das 134 agências e agregá-las (*bottom-up*), ignora-se a influência que os níveis mais agregados — aqui a carteira de crédito da região ou de todo o estado — podem ter na estimação do saldo de crédito de cada agência. De qualquer forma, utiliza-se informação incompleta. Nesse sentido, a reconciliação ótima pode incrementar em grande medida a acurácia das previsões de saldos de crédito de instituições financeiras, permitindo a elaboração do *guidance* e do orçamento de forma mais precisa.

Neste estudo de caso, buscou-se aprimorar as previsões de saldos de empréstimos e financiamentos do Banco do Estado do Espírito Santo através da aplicação da metodologia de combinação ótima via ML. Para o *benchmarking*, foram testadas 7 implementações de métodos de ML (XGB, RF, LGBM, *Elastic Net*, *Lasso*, *ridge* e SVM), 2 métodos analíticos (*bottom-up* e MinT), 3 estratégias de aquisição do conjunto de treinamento para os modelos de ML (*rolling forecast*, *fitted values* e *reduced fitted values*) e 2 estratégias de validação cruzada (aleatória e *rolling origin*).

## Revisão de literatura

A nível macroeconômico, a previsão do agregado de crédito das instituições financeiras é uma preocupação de bancos centrais ao redor do mundo. No Brasil, @bader_modelo_2014 aprimoram o método FAVAR (*Factor Augmented Vector Autoregression*) com uma etapa de análise de correlação canônica para identificar as melhores, em termos de correlação com as variáveis de crédito do SFN, combinações lineares de componentes principais. @colak_tcmb_2019 produzem, a partir de séries filtradas do agregado de crédito, indicadores para monitoramento de períodos de expansão e desaceleração de crédito no setor bancário turco.

Já para níveis abaixo do agregado de crédito, poucos trabalhos foram encontrados. Tangenciando o tema da previsão de saldos de crédito, outros tópicos da economia bancária foram objeto de estudo para previsão de séries temporais. @sezer_financial_2019 produziram revisão de literatura de trabalhos que realizaram previsão de séries temporais financeiras utilizando *deep learning*. @gorodetskaya_machine_2021 propõem o que chamaram de "uma metodologia universal" para aplicação automática de *machine learning* na previsão de séries temporais do setor bancário. Entretanto, nenhum dos trabalhos combinaram estruturas hierárquicas com *machine learning*.

No que diz respeito à previsão de séries temporais em largas hierarquias, @prayoga_top-down_2017 trabalharam na previsão do fluxo de caixa do Banco da Indonésia, utilizando uma hierarquia de 3 níveis. Porém, utilizaram apenas a abordagem *top-down* para reconciliação. Inversamente, @li_hierarchical_2016 compararam dois métodos de *machine learning* para previsão da produção de energia solar no estado da Flórida/EUA, redes neurais artificiais e SVM, em uma abordagem *bottom-up*, caracterizando ambos trabalhos como de nível único.

Aplicações de reconciliação ótima foram realizadas a nível macroeconômico para previsão do PIB australiano, índices de inflação no Reino Unido e México, e previsão de desemprego no Brasil [@athanasopoulos_forecast_2023]. Entretanto, não foram encontradas referências à aplicação de reconciliação ótima para previsão hierárquica de saldos de empréstimos e financiamentos de bancos, tanto a nível individual quanto para o agregado de crédito de uma economia, de forma que este trabalho possa ser o primeiro a aplicar a abordagem de combinação ótima via ML para esse contexto.

## Metodologia {#sec-metodologia_estban}

Os dados usados nesse trabalho são dados terciários obtidos do *datalake* público Base dos Dados [@dahis_data_2022]. A fonte primária são os bancos comerciais e múltiplos com carteira comercial que disponibilizam mensalmente os saldos dos principais verbetes do balancete via documento 4500^[Esses documentos são relatórios eletrônicos obrigatórios demandados pelo Bacen às instituições financeiras que permitem ao regulador o conhecimento minucioso dos bancos e de seus clientes.] ao Banco Central do Brasil, que os compila e publica, agrupados por agência bancária e por município, no relatório ESTBAN — Estatística Bancária Mensal e por Município^[https://www4.bcb.gov.br/fis/cosif/estban.asp?frame=1].

O que compõe os verbetes de crédito, ou seja, os valores das séries temporais a serem trabalhadas, são os saldos de crédito ativo (empréstimos e financiamentos), que correspondem ao principal mais os juros calculados até 59 dias de atraso^[Não são consideradas crédito ativo as operações de crédito liquidadas ou que tenham sido transferidas para prejuízo. São transferidas para prejuízo as operações de crédito em atraso há mais 6 meses após sua classificação de risco em H, que é a mínima [@conselho_monetario_nacional_resolucao_1999].]. Além das estatísticas bancárias, foram obtidas informações de regiões, mesorregiões e microrregiões dos estados, também a partir *datalake* Base dos Dados, com o objetivo de enriquecer a estrutura hierárquica dos dados do ESTBAN, limitada aos municípios.

Em relação à cobertura temporal, utilizou-se dados entre 2003 e 2022. Por brevidade, optou-se por não trabalhar com dados faltantes^[Ver @sec-metodologia_ml.]. Portanto, apesar do Banestes contar com 134 agências, foram mantidas apenas as agências com série completa, ou seja, que já estavam em atividade em 2003.

Por fim, o conjunto de dados conta com séries mensais de saldos de empréstimos e financiamentos de `r with(subset(estban, !is_aggregated(cnpj_agencia)), length(unique(cnpj_agencia)))` agências bancárias, distribuídas por `r with(subset(estban, !is_aggregated(nome)), length(unique(nome)))` municípios, entre `r format(min(estban$ref), "%B de %Y")` e `r format(max(estban$ref), "%B de %Y")`, totalizando `r scales::number(nrow(subset(estban, !is_aggregated(cnpj_agencia) & !is_aggregated(verbete))), big.mark = ".", decimal.mark = ",")` observações. Esses dados então são organizados de forma hierárquica por estado, mesorregião, microrregião, município e agência bancária; e, de forma agrupada, por verbete. Com a estrutura hierárquica e agrupada, o conjunto de dados alcança `r scales::number(nrow(estban))` observações.

```{r}
#| tbl-cap: "Variáveis do banco de dados"
#| label: qdr-variaveis
#| include: false

tibble::tribble(
  ~Variável, ~Descrição,
  "ref", "Data de referência do relatório ESTBAN",
  "nome mesorregiao", "Nome da mesorregião do ES",
  "nome microrregiao", "Nome da microrregião do ES",
  "verbete", "Descrição da rubrica do balancete",
  "nome", "Nome do município",
  "cnpj agencia", "CNPJ da agência bancária",
  "saldo", "Saldo do verbete"
) |>
  kbl(booktabs = TRUE, row.names = FALSE, format = "latex", caption = "Variáveis do banco de dados", label = "qdr-variaveis") |>
  kable_styling(latex_options = c("striped"), font_size = 12) |>
  column_spec(2, width = "10cm")
```

```{r}
#| label: tbl-dataset
#| tbl-cap: "Estrutura do dataset"
#| include: false

estban |>
  tail(13) |>
  transform(
    # switch decimal e big mark
    saldo = scales::number(
      saldo,
      big.mark = ".", decimal.mark = ","
    ),
    # substitui "_" por espaço em todas as colunas
    verbete = gsub("_", " ", verbete),
    municipio = gsub("_", " ", nome),
    microrregiao = gsub("_", " ", nome_microrregiao),
    mesorregiao = gsub("_", " ", nome_mesorregiao),
    agencia = cnpj_agencia
  ) |>
  subset(
    select = c(
      ref,
      mesorregiao,
      microrregiao,
      municipio,
      agencia,
      verbete,
      saldo
    )
  ) |>
  kbl(
    booktabs = TRUE,
    row.names = FALSE
  ) |>
  kable_styling(latex_options = c("striped", "scale_down"))
```

O estado do Espírito Santo está localizado no sudeste brasileiro e é dividido em 78 municípios, que estão agrupados em 4 mesorregiões e 13 microrregiões.

::: {#fig-mapas-es layout-ncol=2}

![Posição no Brasil](img/mapa-brasil.png){#fig-mapa-brasil}

![Micro e mesorregiões](img/mapa-microrregioes.png){#fig-mapa-es}

![Quantidade de agências por município](img/mapa-municipalidades.png){#fig-mapa-municipalidades}

O Estado do Espírito Santo e suas meso e microrregiões.
:::

```{r}
#| label: tbl-microrregioes
#| tbl-cap: "Microrregiões por mesorregião"

aggregate(
  as.character(nome_microrregiao) ~ as.character(nome_mesorregiao),
  data = subset(estban, !is_aggregated(nome_mesorregiao) & !is_aggregated(nome_microrregiao)),
  FUN = function(x) paste(unique(x), collapse = ", ")
) |>
  #substituir "_" por espaço
  transform(
    nome_mesorregiao = tools::toTitleCase(gsub("_", " ", `as.character(nome_mesorregiao)`)),
    nome_microrregiao = tools::toTitleCase(gsub("_", " ", `as.character(nome_microrregiao)`))
  ) |>
  # substituir "Do", "Da" e "De" por "do", "da" e "de"
  transform(
    nome_microrregiao = gsub("Da", "da", gsub("De", "de", gsub("Do", "do", nome_microrregiao)))
  ) |>
  subset(select = c(nome_mesorregiao, nome_microrregiao)) |>
  kbl(
    booktabs = TRUE,
    row.names = FALSE,
    col.names = c(
      "Mesorregião",
      "Microrregiões"
    )
  ) |>
  kable_styling(latex_options = c("striped"))
```

```{r}
#| label: tbl-municipios
#| tbl-cap: "Municípios por microrregião"

aggregate(
  as.character(nome) ~ as.character(nome_microrregiao),
  data = subset(estban, !is_aggregated(nome) & !is_aggregated(nome_microrregiao)),
  FUN = function(x) paste(unique(x), collapse = ", ")
) |>
  # substituir "_" por espaço e adicionar title case
  transform(
    nome_microrregiao = tools::toTitleCase(gsub("_", " ", `as.character(nome_microrregiao)`)),
    nome = tools::toTitleCase(gsub("_", " ", `as.character(nome)`))
  ) |>
  # substituir "Do" e "De" por "do" e "de"
  transform(
    nome = gsub("Da", "da", gsub("De", "de", gsub("Do", "do", nome))),
    nome_microrregiao = gsub("Da", "da", gsub("De", "de", gsub("Do", "do", nome_microrregiao)))
  ) |>
  subset(select = c(nome_microrregiao, nome)) |>
  kbl(
    booktabs = TRUE,
    row.names = FALSE,
    col.names = c(
      "Microrregiões",
      "Municípios"
    )
  ) |>
  kable_styling(latex_options = c("striped"), font_size = 10) |>
  column_spec(2, width = "10cm")
```

Na amostra selecionada para este trabalho, com exceção dos municípios de Colatina e Cachoeiro de Itapemirim, que contam com 2 agências cada, todos os demais municípios no interior contam com uma única agência. Já na microrregião da capital Vitória, se encontram municípios com múltiplas unidades. Considerando que essa microrregião é a mais representativa (@fig-micro), é importante a manutenção do nível hierárquico municipal.

```{r skim}
#| label: tbl-estban
#| tbl-cap: "Contagem de únicos no dataset ESTBAN"
#| include: false

temp = lapply(
  subset(estban, !is_aggregated(cnpj_agencia) & !is_aggregated(verbete))[sapply(estban, function(x) !is.numeric(x))],
  unique
) |>
  lengths()

 # substituir "_" por espaço
names(temp) = gsub("_", " ", names(temp))

temp |>  
  kbl(booktabs = TRUE, col.names = "Únicos") |>
  kable_styling(latex_options = c("striped"))
```

A série temporal do agregado de crédito no Banestes no Espírito Santo é apresentada na @fig-agregado. Em relação à distribuição, a mesorregião Central Espírito-santense concentra mais crédito do que o somatório das demais regiões (@fig-meso), sendo a microrregião de Vitória a responsável por essa concentração (@fig-micro). Na ótica dos verbetes, o crédito para financiamentos é uma pequena fração do saldo de empréstimos, independentemente da mesorregião (@fig-verbetes e @fig-verbete-meso). 

```{r}
#| label: fig-agregado
#| fig-cap: "Série temporal do agregado de crédito do Banestes no ES"

# plot agregado
estban |>
  dplyr::filter(
    is_aggregated(nome),
    is_aggregated(cnpj_agencia),
    is_aggregated(nome_mesorregiao),
    is_aggregated(nome_microrregiao),
    is_aggregated(verbete)
  ) |>
  ggplot(aes(x = ref, y = saldo)) +
  geom_line() +
  scale_y_continuous(labels = scales::number_format(scale = 1 / 1e9)) +
  labs(x = "", y = "Saldo (em R$ bi)") +
  tema
```

```{r}
#| label: fig-meso
#| fig-cap: "Séries temporais do agregado de crédito do Banestes por mesorregião do ES"

# plot agregado
estban |>
  dplyr::filter(
    is_aggregated(nome),
    is_aggregated(cnpj_agencia),
    !is_aggregated(nome_mesorregiao),
    is_aggregated(nome_microrregiao),
    is_aggregated(verbete)
  ) |>
  # substituir "_" por espaço
  transform(
    nome_mesorregiao = tools::toTitleCase(gsub("_", " ", nome_mesorregiao)),
    nome_microrregiao = tools::toTitleCase(gsub("_", " ", nome_microrregiao)),
    nome = tools::toTitleCase(gsub("_", " ", nome)),
    cnpj_agencia = gsub("_", " ", cnpj_agencia),
    verbete = tools::toTitleCase(gsub("_", " ", verbete))
  ) |>
  ggplot(aes(x = ref, y = saldo)) +
  geom_line() +
  scale_y_continuous(labels = scales::number_format(scale = 1 / 1e9)) +
  scale_x_yearmonth(date_labels = "%Y") +
  facet_wrap(~ factor(nome_mesorregiao), scales = "free") +
  labs(x = "", y = "Saldo (em R$ bi)") +
  tema
```

```{r}
#| label: fig-micro
#| fig-cap: "Séries temporais do agregado de crédito do Banestes por microrregião do ES"

# plot agregado
temp = estban |>
  subset(
    is_aggregated(nome)
    & is_aggregated(cnpj_agencia)
    & !is_aggregated(nome_mesorregiao)
    & !is_aggregated(nome_microrregiao)
    & is_aggregated(verbete)
  ) |>
  # substituir "_" por espaço
  transform(
    nome_mesorregiao = gsub("_", " ", nome_mesorregiao),
    nome_microrregiao = gsub("_", " ", nome_microrregiao),
    nome = gsub("_", " ", nome),
    cnpj_agencia = gsub("_", " ", cnpj_agencia),
    verbete = gsub("_", " ", verbete)
  ) |>
  transform(
    nome_mesorregiao = as.character(nome_mesorregiao),
    nome_microrregiao = as.character(nome_microrregiao),
    nome = as.character(nome),
    cnpj_agencia = as.character(cnpj_agencia),
    verbete = as.character(verbete)
  )

plots = temp |>
  by(INDICES = temp$nome_mesorregiao, FUN = function(x) {
    g = ggplot(x, aes(x = ref, y = saldo, color = factor(nome_microrregiao))) +
      geom_line() +
      scale_y_continuous(labels = scales::number_format(scale = 1 / 1e6)) +
      scale_x_yearmonth(date_labels = "%Y") +
      facet_wrap(~ factor(nome_mesorregiao), scales = "free_y") +
      labs(x = "", y = "Saldo (em R$ mi)", color = "microrregião") +
      tema +
      guides(
        color = guide_legend(position = "inside")
      ) +
      theme(
        legend.key.size = unit(0.2, "cm"),
        legend.title = element_text(size = 6),
        legend.text = element_text(size = 6),
        strip.text = element_text(size = 6),
        axis.title.y = element_text(size = 6),
        legend.justification.inside = c(0.08, 0.9)
      )

  })

do.call(gridExtra::grid.arrange, c(plots, ncol = 2))
```

```{r}
#| label: fig-verbetes
#| fig-cap: "Séries temporais dos verbetes no agregado do ES"

# plot por verbetes
estban |>
  dplyr::filter(
    is_aggregated(nome),
    is_aggregated(cnpj_agencia),
    is_aggregated(nome_mesorregiao),
    is_aggregated(nome_microrregiao),
    !is_aggregated(verbete)
  ) |>
  # substituir "_" por espaço
  transform(
    nome_mesorregiao = gsub("_", " ", nome_mesorregiao),
    nome_microrregiao = gsub("_", " ", nome_microrregiao),
    nome = gsub("_", " ", nome),
    cnpj_agencia = gsub("_", " ", cnpj_agencia),
    verbete = gsub("_", " ", verbete)
  ) |>
  ggplot(aes(x = ref, y = saldo, color = factor(verbete))) +
  geom_line() +
  scale_y_continuous(labels = scales::number_format(scale = 1 / 1e9)) +
  labs(x = "", y = "Saldo (em R$ bi)", color = "") +
  tema +
  theme(legend.position = "bottom")
```

```{r}
#| label: fig-verbete-meso
#| fig-cap: "Séries temporais dos verbetes por mesorregião do ES"

# plot por mesorregião
estban |>
  dplyr::filter(
    is_aggregated(nome),
    is_aggregated(cnpj_agencia),
    !is_aggregated(nome_mesorregiao),
    is_aggregated(nome_microrregiao),
    !is_aggregated(verbete)
  ) |>
  # substituir "_" por espaço
  transform(
    nome_mesorregiao = gsub("_", " ", nome_mesorregiao),
    nome_microrregiao = gsub("_", " ", nome_microrregiao),
    nome = gsub("_", " ", nome),
    cnpj_agencia = gsub("_", " ", cnpj_agencia),
    verbete = gsub("_", " ", verbete)
  ) |>
  ggplot(aes(x = ref, y = saldo, color = factor(verbete))) +
  geom_line() +
  facet_wrap(~ factor(nome_mesorregiao), scales = "free_y") +
  scale_y_continuous(labels = scales::number_format(scale = 1 / 1000000)) +
  scale_x_yearmonth(date_labels = "%Y") +
  labs(x = "", y = "Saldo (em R$ mi)", color = "") +
  tema +
  theme(legend.position = "bottom")
```


Em relação às previsões base, uma vez que o foco deste trabalho está no incremento de performance proporcionado pela reconciliação ótima sobre um conjunto de previsões individuais para fora da amostra e não na qualidade destas, elas foram obtidas por meio de métodos simples para previsão de séries temporais, especificamente o algoritmo de Hyndman-Koehler-Snyder-Grose para suavização exponencial (ETS) [@hyndman_state_2002]. Esses métodos são amplamente utilizados na literatura de séries temporais e, portanto, servem como *benchmark* para a avaliação dos métodos de reconciliação ótima.

Os modelos de previsão base foram treinados com dados de 2003 a 2021 e as previsões $h$ passos à frente foram realizadas para o ano de 2022, com $h=\{1, 2, 3, ..., 12\}$ (@fig-modelagem-1). Para verificar a qualidade básica do ajuste, foram realizados testes de Ljung-Box para os resíduos das previsões. Considerando nível de significância de $\alpha=0.05$ e 12 defasagens, cerca de 80% dos modelos não rejeitam a hipótese nula de que os resíduos são descorrelatados. Isso sugere que os modelos de previsão base são, em sua maior parte, adequados para o propósito deste trabalho.

```{r}
#| label: tbl-lb-pvalue
#| tbl-cap: "Resultados do teste de Ljung-box para as previsões base ($\\alpha = 0.05$)"
#| include: true

testes_lb = readRDS("data/estban/previsoes_base/testes_lb.rds") |>
  transform(h_0 = ifelse(lb_pvalue < 0.05, "rejeita", "não rejeita"))

prop.table(table(testes_lb$h_0)) |>
  kbl(booktabs = TRUE, col.names = c("Resultado", "Proporção"), digits = 3) |>
  kable_styling(latex_options = c("striped"))
```

O *software* usado foi o R [@R-base]. As previsões base, as previsões hierárquicas de nível único e a reconciliação ótima através da MinT foram realizados com os pacotes {fable} [@R-fable] e {fabletools} [@R-fabletools]. Já metodologia de *machine learning* (reamostragem, otimização de hiperparâmetros, treino e predição) foi executada com o pacote {mlr3} [@R-mlr3] e suas extensões. As implementações do XGB e LGBM utilizadas foram as de seus pacotes homônimos {xgboost} [@R-xgboost] e {ligthgbm} [@R-lightgbm]. Para as configurações de *Elastic Net* foi utilizado o pacote {glmnet} [@R-glmnet]. Já para os métodos de RF e SVM foram usados os pacotes {ranger} [@R-ranger] e {e1071} [@meyer__aut_e1071_2023], respectivamente.

Boa parte dos métodos de ML são altamente parametrizáveis, sendo sua performance de generalização sensível à escolha de seus hiperparâmetros. Para os métodos avaliados neste trabalho, os conjuntos de hiperparâmetros a serem otimizados e seus respectivos espaços de busca são apresentados no Apêndice \ref{apendice_hiperparametros}, de acordo com a recomendação em @bischl_hyperparameter_2021.

Para a otimização, foram utilizadas dois calibradores: (i) busca em grade, com resolução de 10 combinações, mais custoso em tempo de processamento, para os métodos de menor quantidade de hiperparâmetros a serem otimizados; e (ii) otimização bayesiana, na configuração padrão do pacote {mlr3MBO} [@R-mlr3mbo], mais eficiente para os métodos com maior quantidade de hiperparâmetros. As estratégias de reamostragem utilizadas foram a validação cruzada *k-fold* aleatória e a *k-fold* com origem móvel, ambas com $k=10$. Para a segunda, utilizou-se janela mínima de 5 anos e $|\Omega^b| = 12$, compatível com o horizonte de previsão para o experimento.

A otimização bayesiana foi usada em todos os métodos, exceto no *elastic net*, uma vez que apenas um (no caso do *Lasso* e *ridge*) ou dois hiperparâmetros foram otimizados. A medida de performance utilizada para a otimização foi a raiz do erro quadrático médio (*root mean squared error* — RMSE).

Por fim, cada modelo foi calibrado individualmente, ou seja, cada série $y_m$ possui um conjunto de hiperparâmetros otimizados para cada um dos 7 métodos de ML empregados para combinação ótima.

<!-- TODO: gráficos de previsão para fora da amostra -->

```{r}
#| label: fig-previsoes
#| fig-cap: Previsões para fora da amostra
#| eval: false

# caminho
caminho = list.files("data/estban/preds_ml/preds/rolling_forecast", full.names = TRUE) |>
  (function(arquivo) subset(x = arquivo, !grepl("resumo", x = arquivo) & grepl("preds", x = arquivo)))()

# extrair do caminho o que está entre "preds_" e ".RDS"
nomes = gsub(".*preds_(.*?)\\.RDS.*", "\\1", caminho)

# juntar previsões
predicoes = lapply(readRDS(caminho[1])[[1]], function(x) {
  data.table::as.data.table(x)[[3]]
}) |>
  stack() |>
  transform(ref = 1:12)

colnames(predicoes) = c("preds", "nomes", "ref")
```

## Resultados

As tabelas a seguir apresentam os resultados obtidos para o experimento. As métricas utilizadas foram a Raiz do Erro Médio Escalado Quadrático (*Root Mean Squared Scaled Error* — RMSSE) e o erro médio escalado absoluto (*Mean Absolute Scaled Error* — MASE), que permitem a comparação independentemente da escala — essencial no contexto de previsão hierárquica, onde a agregação torna natural a diferença de escala.

A @tbl-estban-results-analiticos contém as medidas de acurácia RMSSE e MASE para os métodos analíticos de reconciliação ótima BU e MinT, e para as previsões base, ou seja, sem aplicar qualquer método de reconciliação. A primeira coluna especifica o método utilizado, enquanto as demais colunas apresentam a média da performance em cada nível de agregação.

As @tbl-estban-results-ml-rolling, @tbl-estban-results-ml-fitted e @tbl-estban-results-ml-refit reportam as medidas de acurácia para os métodos de reconciliação ótima baseados em ML utilizando validação cruzada aleatória, enquanto as @tbl-estban-results-ml-rolling-fcv, @tbl-estban-results-ml-fitted-fcv e @tbl-estban-results-ml-refit-fcv exibem os resultado para validação cruzada com origem móvel. Já as @tbl-estban-tempo-ml e @tbl-estban-tempo-ml-fcv reportam o tempo de processamento para as etapas de calibragem, treino e predição desses métodos^[Os métodos analíticos não tiveram seu tempo de processamento medidos porque executam quase que instantaneamente, já sinalizando uma vantagem para esses métodos.].

Em geral, os métodos baseados em árvore, além de requererem maior tempo de processamento devido a sua complexidade no espaço de hiperparâmetros, também tenderam a perder qualidade conforme suas previsões são agregadas para formação dos níveis superiores da hierarquia. Contrariamente, os métodos de regressão regularizada se mostraram mais robustos à agregação.

Nas tabelas a seguir, **negrito** indica a melhor performance entre os métodos para aquele determinado nível de agregação, e \underline{sublinhado} indica que aquele método de ML superou o método analítico de melhor performance naquele nível de agregação.

Nos resultados para a base de dados Estban estão incluídas as médias de performance para cada nível hierárquico e agrupado. As colunas "agregado", "mesorregiao", "microrregiao", "municipio" e "agencia", fazem referência à estrutura hierárquica, ou seja, tratam o verbete de forma agregada. Já as colunas "verbete", "bottom" e "hierarquia", incluem também a estrutura agrupada, tratando o verbete de forma desagregada. Detalhadamente:

- Agregado: performance do método para a série que representa o total, com os verbetes agregados, (@fig-agregado).
- Mesorregião: a média das performances do método para as séries do agregado de cada mesorregião, com os verbetes agregados (@fig-meso).
- Microrregião: a média das performances do método para as séries do agregado de cada microrregião, com os verbetes agregados (@fig-micro).
- Município: a média das performances do método para as séries do agregado de cada município, com os verbetes agregados.
- Agência: a média das performances do método para as séries de cada agência, com os verbetes agregados.
- Verbete: a média das performances do método para as séries de cada verbete, para o total da hierarquia (@fig-verbetes).
- Bottom: a média das performances do método para as séries do nível mais desagregado, ou seja, verbete por agência.
- Hierarquia: a média das performances do método para todas as séries, agregadas e desagregadas.

```{r}
#| tbl-cap: "Resultados Estban: Acurácia dos métodos analíticos de reconciliação"
#| label: tbl-estban-results-analiticos

rbind(
  readRDS("data/estban/preds_analitico/acuracia_analiticos.rds")[["rmsse"]],
  readRDS("data/estban/preds_analitico/acuracia_analiticos.rds")[["mase"]]
 ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 3) |>
  pack_rows("MASE", 4, 6) |>
  kable_styling(latex_options = c("striped", "scale_down")) |>
  # município
  column_spec(
    5,
    bold = c(T, F, F, T, F, F)
  ) |>
  # agência
  column_spec(
    6,
    bold = c(F, F, T, T, F, F)
  ) |>
  # verbete
  column_spec(
    7,
    bold = c(F, F, T, F, F, T)
  ) |>
  # bottom
  column_spec(
    8,
    bold = c(F, F, T, F, T, F)
  ) |>
  # hierarquia
  column_spec(
    9,
    bold = c(F, F, T, F, F, T)
  )
```

```{r}
#| tbl-cap: "Resultados Estban: Acurácia dos métodos de ML de reconciliação, estratégia *rolling forecast* + validação cruzada aleatória"
#| label: tbl-estban-results-ml-rolling

rbind(
  readRDS("data/estban/preds_ml/preds/rolling_forecast/resumo.RDS")[["rmsse"]],
  readRDS("data/estban/preds_ml/preds/rolling_forecast/resumo.RDS")[["mase"]]
 ) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped", "scale_down")) |>
  # agregado
  column_spec(
    2,
    underline = c(F, F, F, F, F, F, F, F, T, F, F, F, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados Estban: Acurácia dos métodos de ML de reconciliação, estratégia *rolling forecast* + validação cruzada com origem móvel"
#| label: tbl-estban-results-ml-rolling-fcv

rbind(
  readRDS("data/estban/preds_ml/preds/rolling_forecast/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/estban/preds_ml/preds/rolling_forecast/resumo_fcv.RDS")[["mase"]]
 ) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped", "scale_down")) |>
  # agregado
  column_spec(
    2,
    underline = c(F, F, F, F, F, F, F, F, T, F, F, F, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados Estban: Acurácia dos métodos de ML de reconciliação, estratégia *fitted values* + validação cruzada aleatória"
#| label: tbl-estban-results-ml-fitted

rbind(
  readRDS("data/estban/preds_ml/preds/fitted_base/resumo.RDS")[["rmsse"]],
  readRDS("data/estban/preds_ml/preds/fitted_base/resumo.RDS")[["mase"]]
 ) |>
 # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped", "scale_down"))
```

```{r}
#| tbl-cap: "Resultados Estban: Acurácia dos métodos de ML de reconciliação, estratégia *fitted values* + validação cruzada com origem móvel"
#| label: tbl-estban-results-ml-fitted-fcv

rbind(
  readRDS("data/estban/preds_ml/preds/fitted_base/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/estban/preds_ml/preds/fitted_base/resumo_fcv.RDS")[["mase"]]
 ) |>
 # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped", "scale_down"))
```

```{r}
#| tbl-cap: "Resultados Estban: Acurácia dos métodos de ML de reconciliação, estratégia *reduced fitted values* + validação cruzada aleatória"
#| label: tbl-estban-results-ml-refit

rbind(
  readRDS("data/estban/preds_ml/preds/one-step-ahead/resumo.RDS")[["rmsse"]],
  readRDS("data/estban/preds_ml/preds/one-step-ahead/resumo.RDS")[["mase"]]
 ) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped", "scale_down")) |>
  # agregado
  column_spec(
    2,
    bold = c(T, F, F, F, F, F, F, T, F, F, F, F, F, F),
    underline = c(F, T, F, F, T, F, F, F, T, F, F, T, F, F)
  ) |>
  # mesorregião
  column_spec(
    3,
    underline = c(F, T, F, F, F, F, F, F, T, F, F, F, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados Estban: Acurácia dos métodos de ML de reconciliação, estratégia *reduced fitted values* + validação cruzada com origem móvel"
#| label: tbl-estban-results-ml-refit-fcv

rbind(
  readRDS("data/estban/preds_ml/preds/one-step-ahead/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/estban/preds_ml/preds/one-step-ahead/resumo_fcv.RDS")[["mase"]]
 ) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped", "scale_down")) |>
  # agregado
  column_spec(
    2,
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F),
    underline = c(T, F, F, F, T, F, F, T, F, F, F, T, F, F)
  ) |>
  # mesorregião
  column_spec(
    3,
    bold = c(F, T, F, F, F, F, F, F, T, F, F, F, F, F)
  ) |>
  # microrregião
  column_spec(
    4,
    bold = c(F, T, F, F, F, F, F, F, T, F, F, F, F, F),
    underline = c(T, F, F, F, F, F, F, T, F, F, F, F, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados Estban: Tempo de processamento dos métodos de ML (em horas). Validação cruzada aleatória"
#| label: tbl-estban-tempo-ml

# obtendo tempo de processamento
tempo = lapply(c("one-step-ahead", "fitted_base", "rolling_forecast"), function(tipo) {
  lapply(c("xgb", "ranger", "glmnet", "lasso", "ridge", "svm", "lightgbm"), function(learner) {
    preds = readRDS(paste0("data/estban/preds_ml/preds/", tipo, "/preds_", learner, ".RDS"))[[2]]
    return(preds)
  })
})

# juntando predições em um único dataframe
tempo = do.call(rbind, tempo)

# convertendo colunas para numérico
tempo = sapply(as.data.frame(tempo), function(x) as.numeric(x))

# nomeando df
colnames(tempo) = c("xgb", "random forest", "elastic net", "lasso", "ridge", "svm", "lightgbm")
rownames(tempo) = c("reduced fitted values", "fitted values", "rolling forecast")

# tabela
tempo |>
  kbl(booktabs = TRUE, digits = 3) |>
  kable_styling(latex_options = c("striped"))
```

```{r}
#| tbl-cap: "Resultados Estban: Tempo de processamento dos métodos de ML (em horas). Validação cruzada com origem móvel"
#| label: tbl-estban-tempo-ml-fcv

# obtendo tempo de processamento
tempo = lapply(c("one-step-ahead", "fitted_base", "rolling_forecast"), function(tipo) {
  lapply(c("xgb", "ranger", "glmnet", "lasso", "ridge", "svm", "lightgbm"), function(learner) {
    preds = readRDS(paste0("data/estban/preds_ml/preds/", tipo, "/preds_", learner, "_fcv.RDS"))[[2]]
    return(preds)
  })
})

# juntando predições em um único dataframe
tempo = do.call(rbind, tempo)

# convertendo colunas para numérico
tempo = sapply(as.data.frame(tempo), function(x) as.numeric(x))

# nomeando df
colnames(tempo) = c("xgb", "random forest", "elastic net", "lasso", "ridge", "svm", "lightgbm")
rownames(tempo) = c("reduced fitted values", "fitted values", "rolling forecast")

# tabela
tempo |>
  kbl(booktabs = TRUE, digits = 3) |>
  kable_styling(latex_options = c("striped"))
```

Para a base de dados Estban, não houve uma combinação de método e estratégia que fosse consistentemente melhor ao longo de todos os níveis de agregação. Portanto, a escolha do método e da estratégia a serem utilizados dependerá do objetivo do pesquisador^[Se o objetivo é a elaboração de *guidance*, por exemplo, o pesquisador deve preferir o método e estratégia que geram as  previsões mais precisas para o nível agregado. Já para elaboração de metas individuais, os níveis individuais ou regionais podem ser preferíveis.].

Para os níveis ao topo da hierarquia, os métodos de ML se mostraram a melhor opção para estimação. No nível agregado, o *elastic net* na estratégia *reduced fitted values* + validação cruzada aleatória (@tbl-estban-results-ml-refit) se mostrou a melhor opção para a estimação do agregado, com 91% de ganho de performance sobre o MinT (@tbl-estban-results-analiticos), em termos de RMSSE, ou seja, o MinT obteve quase o dobro do erro do *elastic net*. Para os níveis de meso e microrregião, o *elastic net* na configuração *Lasso* e na estratégia *reduced fitted values* + validação cruzada com origem móvel obteve as melhores performances, sendo 10% melhor que BU para mesorregião e 6% melhor que o BU para microrregião (@tbl-estban-results-ml-refit-fcv). Por outro lado, tanto nos níveis hierárquicos abaixo, quanto no nível agrupado, os métodos de ML não foram capazes de superar os métodos analíticos.

Os resultados se mostraram bastante sensíveis à estratégia utilizada. Na métrica RMSS, nenhum método utilizando as estratégias *rolling forecast* e *fitted values* foi capaz de superar o MinT em qualquer nível de agregação, enquanto na estratégia *reduced fitted values* os métodos de ML mostraram ganhos de performance nas três configurações do *elastic net* em ao menos um nível de agregação.

Em geral, as medidas RMSSE e MASE se mostraram bastante correlacionadas, com os métodos que obtiveram melhor performance em uma métrica também obtendo melhor performance na outra. A exceção foi para o método *lasso* na estratégia *rolling forecast*, que obteve performance melhor que o MinT, em termos de MASE, para o nível agregado, sendo o único resultado positivo para a estratégia *rolling forecast*.

## Conclusão

Neste estudo de caso, foram apresentados experimentos de reconciliação ótima para séries temporais hierárquicas e agrupadas, utilizando métodos analíticos e de ML com o objetivo de obter coerência e ganhos de acurácia nas previsões de saldos de empréstimos e financiamentos do Banco do Estado do Espírito Santo. Pesquisas anteriores já haviam mostrado que a reconciliação ótima pode trazer ganhos de acurácia, e que métodos de ML podem ser competitivos em relação aos métodos analíticos.

Este trabalho trouxe, além dos métodos de RF e XGB já trabalhados em @spiliotis_hierarchical_2021, o método de regressão regularizada *elastic net*, também em suas configurações *Lasso* e *ridge*, e o SVM, além de avaliar outra implementação de *gradient boosting*, o LGBM. Paralelamente, este trabalho propôs duas estratégias alternativas para aquisição do conjunto de treinamento dos modelos de ML e também avaliou a estratégia de validação cruzada com origem móvel para a tarefa de combinação ótima.

Os resultados obtidos para a base de dados Estban mostraram, primeiramente, que não houve uma combinação de método e estratégia que obtivesse melhor performance de maneira consistente ao longo de todos os níveis hierárquicos. Dessa forma, a escolha do método e da estratégia a serem utilizados dependerá do objetivo do pesquisador. Para os níveis ao topo da hierarquia, a combinação correta de método e estratégia de estimação (*elastic net* + *reduced fitted values*) gerou 91% de ganho de performance no nível mais agregado, o que imputaria a instituição financeira maior precisão para tomada de decisão e planejamento estratégico, além de sinalizar maior confiança nas estimativas comunicadas ao mercado e aos investidores.

Por outro lado, para os níveis inferiores na estrutura hierárquica (municípios e agências), os métodos analíticos se mostraram a melhor escolha. Isto sugere que, com o objetivo de elaboração de metas individuais — seja para as agências ou superintendências regionais —, os métodos analíticos ainda são preferíveis, enquanto para microrregiões, mesorregiões e total, os métodos de ML demonstram maior acurácia.

Os restultados também mostram que o resultado da reconciliação ótima é sensível à estratégia de aquisição do conjunto de treinamento dos modelos de ML. Apenas na estratégia *reduced fitted values*, proposta neste trabalho, que os ganhos de performance foram observados, enquanto nas estratégias *rolling forecast* e *fitted values* os métodos de ML não foram capazes de superar os métodos analíticos. Em relação às estratégias de reamostragem, a validação cruzada com origem móvel se mostrou eficiente para otimização de hiperparâmetros na tarefa de combinação ótima, tendo gerado os modelos mais precisos para os níveis de meso e microrregião.

Contrariamente aos resultados de @spiliotis_hierarchical_2021, as implementações de ML baseados em árvore de decisão (XGB, RF e LGBM) não foram capazes de superar os métodos analíticos em nenhum nível de agregação.

# PREVISÕES PARA O TURISMO DOMÉSTICO AUSTRALIANO {#sec-tourism}

## Metodologia

Com o objetivo de avaliar a metodologia proposta em diferentes contextos e compará-la aos trabalhos anteriores, foi realizado um segundo estudo de caso com as bases de dados *overnight trips* e *domestic visitor nights*, ambas relacionadas ao setor de turismo australiano. 

A base de dados *overnight trips* consiste na quantidade trimestral de viagens domésticas que duram uma noite na Austrália entre 1998 e 2017. A estrutura é hierárquica e agrupada, composta por 3 níveis hierárquicos --- \emph{State} (Estados), \emph{Region} (Regiões) e total ---, e agrupado por \emph{Purpose}
(Propósito). Essa base de dados foi utilizada em @athanasopoulos_hierarchical_2009 para comparar a performance da combinação ótima contra BU e algumas variações de TD e é também utilizada em @hyndman_forecasting_2021 como recurso didático no capítulo de previsão hierárquica, onde uma análise exploratória dos dados pode ser encontrada. Esta base de dados está disponível no pacote {tsibble} [@R-tsibble].

Para esse estudo de caso, o *benchmark* será composto pelos métodos analíticos de reconciliação BU e MinT, e pelos métodos de ML *elastic net*, *Lasso*, *ridge*, RF, XGB, LGBM e SVM. Serão avaliadas a estratégia *rolling forecast* e as duas propostas neste trabalho, *fitted values* e *reduced fitted values*, com validação cruzada aleatória e com origem móvel. 

Os modelos de previsões base, obtidos pelo algoritmo de Hyndman-Koehler-Snyder-Grose para ETS, foram treinados com dados entre 1998 a 2015 e as previsões $h$ passos à frente forma realizadas para os anos de 2016 e 2017, com $h = \{1, 2, \ldots, 8\}$. Foram realizados testes de Ljung-Box para verificar a presença de autocorrelação nos resíduos das previsões base, com $\alpha = 0.05$. Considerando que 99% das séries não apresentam autocorrelação serial nos resíduos, espera-se que a validação cruzada aleatória seja suficiente para avaliação dos modelos [@bergmeir_note_2018].

```{r}
#| label: tbl-lb-pvalue-tourism
#| tbl-cap: "Resultados do teste de Ljung-box para as previsões base ($\\alpha = 0.05$) da base de dados *overnight trips*"
#| include: true

testes_lb = readRDS("data/tourism/previsoes_base/testes_lb.rds") |>
  transform(h_0 = ifelse(lb_pvalue < 0.05, "rejeita", "não rejeita"))

prop.table(table(testes_lb$h_0)) |>
  kbl(booktabs = TRUE, col.names = c("Resultado", "Proporção"), digits = 3) |>
  kable_styling(latex_options = c("striped"))
```

A base de dados *domestic visitor nights* registra o número total de noites que os visitantes domésticos passam em viagens na Austrália. A estrutura desse dataset é estritamente hierárquica, composta por 4 níveis hierárquicos: \emph{State} (Estados), \emph{Zones} (Zonas), \emph{Region} (Regiões) e total. Essa base de dados foi utilizada em trabalhos como @wickramasuriya_optimal_2019, @kourentzes_cross-temporal_2019 e @spiliotis_hierarchical_2021. Por brevidade, dispensa-se a análise exploratória dos dados, que pode ser encontrada nos artigos supracitados.

Por essa base de dados ser estritamente hierárquica, além dos métodos empregados na base *overnight trips*, também será avaliado o método TD por proporções de previsão [@athanasopoulos_hierarchical_2009], disponível no pacote {hts} [@R-hts]. Todas as estratégias supracitadas para a base *overnight trips* serão avaliadas também para essa base de dados.

Os modelos de previsões base foram treinados com dados entre 1998 a 2015 e as previsões $h$ passos à frente foram realizadas para os anos de 2016 e 2017, com $h = \{1, 2, \ldots, 12\}$. Para manter alguma semelhança com o experimento de @spiliotis_hierarchical_2021, os modelos foram treinados utilizando o algoritmo de Hyndman-Khandakar para modelagem Arima automática [@hyndman_automatic_2008]. Assim como na base *overnight trips*, foram realizados testes de Ljung-Box para verificar a presença de autocorrelação nos resíduos das previsões base, com $\alpha = 0.05$. Com cerca de 14% dos modelos apresentando autocorrelação serial nos resíduos, espera-se algum ganho de performance com a validação cruzada com origem móvel.

```{r}
#| label: tbl-lb-pvalue-tourism-monthly
#| tbl-cap: "Resultados do teste de Ljung-box para as previsões base ($\\alpha = 0.05$) da base de dados *visitor nights*"
#| include: true

testes_lb = readRDS("data/tourism_monthly/previsoes_base/testes_lb.rds") |>
  transform(h_0 = ifelse(lb_pvalue < 0.05, "rejeita", "não rejeita"))

prop.table(table(testes_lb$h_0)) |>
  kbl(booktabs = TRUE, col.names = c("Resultado", "Proporção")) |>
  kable_styling(latex_options = c("striped"))
```

Para os métodos de ML, foram otimizados o conjunto de hiperparâmetros em seus respectivos espaços de busca informados no Apêndice \ref{apendice_hiperparametros}. Cada série $y_m$ possui um conjunto de hiperparâmetros otimizados para cada um dos 7 métodos de ML empregados. Para a base de dados *overnight trips*, de dados trimestrais, a estratégia de reamostragem por validação cruzada com origem móvel contou com $|\Omega^a| \geq 12$ e $|\Omega^b| = 4$, ou seja, janela mínima de 12 observações (3 anos) e horizonte de 4 observações (1 ano). Já para a base de dados *domestic visitor nights*, de dados mensais, a estratégia de reamostragem contou com $|\Omega^a| \geq 60$ e $|\Omega^b|=12$. Demais pontos em relação a *softwares* e configurações utilizadas seguem os mesmos padrões da seção @sec-metodologia_estban.

## Resultados e discussão

Primeiramente, para a base de dados *overnight trips*, a tabela @tbl-tourism-results-analiticos contém as medidas RMSSE e MASE para os métodos analíticos de reconciliação, BU e MinT, e para as previsões base, ou seja, sem aplicar qualquer método de reconciliação. A primeira coluna especifica o método utilizado, enquanto as demais colunas apresentam a média da performance em cada nível de agregação.

As @tbl-tourism-results-ml-rolling, @tbl-tourism-results-ml-fitted e @tbl-tourism-results-ml-refit reportam as medidas de acurácia para os métodos de reconciliação ótima baseados em ML utilizando validação cruzada aleatória, enquanto as @tbl-tourism-results-ml-rolling-fcv, @tbl-tourism-results-ml-fitted-fcv e @tbl-tourism-results-ml-refit-fcv reportam as medidas de acurácia para os métodos de reconciliação ótima baseados em ML utilizando validação cruzada com origem móvel. As @tbl-tourism-tempo-ml e @tbl-tourism-tempo-ml-fcv apresentam o tempo de processamento dos métodos de ML em horas, para validação cruzada aleatória e com origem móvel, respectivamente.

Para essa base de dados, a aplicação da reconciliação foi capaz de melhorar a qualidade da previsão em relação às previsões base em todos os níveis de agregação, com os métodos de ML alcançando as melhores performances em todos os níveis de agregação, exceto no mais desagregado, onde o MinT obteve a melhor performance. O SVM, na estratégia *rolling forecast* + validação cruzada aleatória, se mostrou a melhor combinação dentre os métodos de ML, alcançando 93% de incremento de performance em relação ao MinT, em termos de RMSSE, ou seja, o MinT teve quase o dobro do erro que o SVM.

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Acurácia dos métodos analíticos de reconciliação"
#| label: tbl-tourism-results-analiticos

rbind(
  readRDS("data/tourism/preds_analitico/acuracia_analiticos.rds")[["rmsse"]],
  readRDS("data/tourism/preds_analitico/acuracia_analiticos.rds")[["mase"]]
 ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 3) |>
  pack_rows("MASE", 4, 6) |>
  kable_styling(latex_options = c("striped")) |>
  # bottom
  column_spec(
    6,
    bold = c(F, F, T, F, F, T)
  ) |>
  # hierarquia
  column_spec(
    7,
    bold = c(F, F, T, F, F, T)
  )
```

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Acurácia dos métodos de ML de reconciliação, estratégia *rolling forecast* + validação cruzada aleatória"
#| label: tbl-tourism-results-ml-rolling

rbind(
  readRDS("data/tourism/preds_ml/preds/rolling_forecast/resumo.RDS")[["rmsse"]],
  readRDS("data/tourism/preds_ml/preds/rolling_forecast/resumo.RDS")[["mase"]]
 ) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped")) |>
  # agregado
  column_spec(
    2,
    bold = c(F, F, F, F, F, T, F, F, F, F, F, F, T, F),
    underline = c(F, F, F, F, T, F, F, F, F, F, F, T, F, F)
  ) |>
  # state
  column_spec(
    3,
    bold = c(F, F, F, F, F, T, F, F, F, F, F, F, T, F),
    underline = c(F, F, F, F, T, F, F, F, F, F, F, T, F, F)
  ) |>
  # purpose
  column_spec(
    5,
    bold = c(F, F, F, F, F, T, F, F, F, F, F, F, T, F),
    underline = c(F, F, F, F, T, F, F, F, F, F, F, T, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Acurácia dos métodos de ML de reconciliação, estratégia *rolling forecast* + validação cruzada com origem móvel"
#| label: tbl-tourism-results-ml-rolling-fcv

rbind(
  readRDS("data/tourism/preds_ml/preds/rolling_forecast/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/tourism/preds_ml/preds/rolling_forecast/resumo_fcv.RDS")[["mase"]]
 ) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped")) |>
  # agregado
  column_spec(
    2,
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F),
    underline = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  ) |>
  # state
  column_spec(
    3,
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F),
    underline = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  ) |>
  # purpose
  column_spec(
    5,
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F),
    underline = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Acurácia dos métodos de ML de reconciliação, estratégia *fitted values* + validação cruzada aleatória"
#| label: tbl-tourism-results-ml-fitted

rbind(
  readRDS("data/tourism/preds_ml/preds/fitted_base/resumo.RDS")[["rmsse"]],
  readRDS("data/tourism/preds_ml/preds/fitted_base/resumo.RDS")[["mase"]]
 ) |>
 # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped")) |>
  # state
  column_spec(
    3,
    underline = c(F, F, F, F, T, T, F, F, F, F, F, T, T, F)
  ) |>
  # region
  column_spec(
    4,
    bold = c(F, F, F, F, T, F, F, F, F, F, F, T, F, F)
  ) |>
  # purpose
  column_spec(
    5,
    underline = c(F, F, F, F, F, F, F, F, F, F, F, T, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Acurácia dos métodos de ML de reconciliação, estratégia *fitted values* + validação cruzada com origem móvel"
#| label: tbl-tourism-results-ml-fitted-fcv

rbind(
  readRDS("data/tourism/preds_ml/preds/fitted_base/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/tourism/preds_ml/preds/fitted_base/resumo_fcv.RDS")[["mase"]]
 ) |>
 # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped")) |>
  # state
  column_spec(
    3,
    underline = c(F, F, F, F, T, F, T, F, F, F, F, T, F, T)
  ) |>
  # region
  column_spec(
    4,
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  ) |>
  # purpose
  column_spec(
    5,
    underline = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Acurácia dos métodos de ML de reconciliação, estratégia *reduced fitted values* + validação cruzada aleatória"
#| label: tbl-tourism-results-ml-refit

rbind(
  readRDS("data/tourism/preds_ml/preds/one-step-ahead/resumo.RDS")[["rmsse"]],
  readRDS("data/tourism/preds_ml/preds/one-step-ahead/resumo.RDS")[["mase"]]
 ) |>
  # aplicando notação científica quando necessário
  sapply(function(x) prettyNum(x, digits = 4)) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped"))
```

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Acurácia dos métodos de ML de reconciliação, estratégia *reduced fitted values* + validação cruzada com origem móvel"
#| label: tbl-tourism-results-ml-refit-fcv

rbind(
  readRDS("data/tourism/preds_ml/preds/one-step-ahead/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/tourism/preds_ml/preds/one-step-ahead/resumo_fcv.RDS")[["mase"]]
 ) |>
  # aplicando notação científica quando necessário
  sapply(function(x) prettyNum(x, digits = 4)) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped"))
```

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Tempo de processamento dos métodos de ML (em horas). Validação cruzada aleatória"
#| label: tbl-tourism-tempo-ml

# obtendo tempo de processamento
tempo = lapply(c("fitted_base", "rolling_forecast", "one-step-ahead"), function(tipo) {
  lapply(c("xgb", "ranger", "glmnet", "lasso", "ridge", "svm", "lightgbm"), function(learner) {
    preds = readRDS(paste0("data/tourism/preds_ml/preds/", tipo, "/preds_", learner, ".RDS"))[[2]]
    return(preds)
  })
})

# juntando predições em um único dataframe
tempo = do.call(rbind, tempo)

# convertendo colunas para numérico
tempo = sapply(as.data.frame(tempo), function(x) as.numeric(x))

# nomeando df
colnames(tempo) = c("xgb", "random forest", "elastic net", "lasso", "ridge", "svm", "lightgbm")
rownames(tempo) = c("fitted values", "rolling forecast", "reduced fitted values")

# tabela
tempo |>
  kbl(booktabs = TRUE, digits = 3) |>
  kable_styling(latex_options = c("striped"))
```

```{r}
#| tbl-cap: "Resultados para *overnight trips*: Tempo de processamento dos métodos de ML (em horas). Validação cruzada com origem móvel"
#| label: tbl-tourism-tempo-ml-fcv

# obtendo tempo de processamento
tempo = lapply(c("fitted_base", "rolling_forecast", "one-step-ahead"), function(tipo) {
  lapply(c("xgb", "ranger", "glmnet", "lasso", "ridge", "svm", "lightgbm"), function(learner) {
    preds = readRDS(paste0("data/tourism/preds_ml/preds/", tipo, "/preds_", learner, "_fcv.RDS"))[[2]]
    return(preds)
  })
})

# juntando predições em um único dataframe
tempo = do.call(rbind, tempo)

# convertendo colunas para numérico
tempo = sapply(as.data.frame(tempo), function(x) as.numeric(x))

# nomeando df
colnames(tempo) = c("xgb", "random forest", "elastic net", "lasso", "ridge", "svm", "lightgbm")
rownames(tempo) = c("fitted values", "rolling forecast", "reduced fitted values")

# tabela
tempo |>
  kbl(booktabs = TRUE, digits = 3) |>
  kable_styling(latex_options = c("striped"))
```

Em relação à base de dados *visitors night*, a tabela @tbl-tourism-monthly-results-analiticos contém as medidas RMSSE e MASE para os métodos analíticos de reconciliação, BU, TD e MinT, e para as previsões base, ou seja, sem aplicar qualquer método de reconciliação. Como nas demais tabelas, a primeira coluna especifica o método utilizado, enquanto as outras apresentam a média da performance em cada nível de agregação.

As @tbl-tourism-monthly-results-ml-rolling, @tbl-tourism-monthly-results-ml-fitted e @tbl-tourism-monthly-results-ml-refit reportam as medidas de acurácia para os métodos de reconciliação ótima baseados em ML utilizando validação cruzada aleatória, enquanto as @tbl-tourism-monthly-results-ml-rolling-fcv, @tbl-tourism-monthly-results-ml-fitted-fcv e @tbl-tourism-monthly-results-ml-refit-fcv reportam as medidas de acurácia para os métodos de reconciliação ótima baseados em ML utilizando validação cruzada com origem móvel. As @tbl-tourism-monthly-tempo-ml e @tbl-tourism-monthly-tempo-ml-fcv apresentam o tempo de processamento dos métodos de ML em horas, para validação cruzada aleatória e com origem móvel, respectivamente.

Para essa base de dados, os métodos de combinação ótima obtiveram a melhor performance em 2 dos 4 níveis hierárquicos, os mais desagregados (*zone* e *region*), enquanto o *top-down* por proporções de previsão foi o método mais preciso nos níveis agregado e *state*. Dentre os métodos de ML, o *elastic net* em sua configuração *ridge*, na estratégia *fitted values* + validação cruzada aleatória, se mostrou a melhor combinação, com 3% de incremento de performance em relação ao MinT no nível *zone* e 2% no nível *region*, também alcanando a melhor média no indicador geral da hierarquia.

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Acurácia dos métodos analíticos de reconciliação"
#| label: tbl-tourism-monthly-results-analiticos

td_data = readRDS("data/tourism_monthly/preds_analitico/acuracia_top_down.rds")

td_rmsse = td_data |> subset(rownames(td_data) == "RMSSE")

td_mase = td_data |> subset(rownames(td_data) == "MASE")

rbind(
  readRDS("data/tourism_monthly/preds_analitico/acuracia_analiticos.rds")[["rmsse"]],
  td_rmsse,
  readRDS("data/tourism_monthly/preds_analitico/acuracia_analiticos.rds")[["mase"]],
  td_mase
 ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 4) |>
  pack_rows("MASE", 5, 8) |>
  kable_styling(latex_options = c("striped")) |>
  # agregado
  column_spec(
    2,
    bold = c(F, F, F, T, F, F, F, T)
  ) |>
  # state
  column_spec(
    3,
    bold = c(F, F, F, T, F, F, F, T)
  )
```

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Acurácia dos métodos de ML de reconciliação, estratégia *rolling forecast* + validação cruzada aleatória"
#| label: tbl-tourism-monthly-results-ml-rolling

rbind(
  readRDS("data/tourism_monthly/preds_ml/preds/rolling_forecast/resumo.RDS")[["rmsse"]],
  readRDS("data/tourism_monthly/preds_ml/preds/rolling_forecast/resumo.RDS")[["mase"]]
 ) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped")) |>
  # zone
  column_spec(
    4,
    underline = c(F, F, F, T, F, F, F, F, F, F, T, F, F, F)
  ) |>
  # region
  column_spec(
    5,
    underline = c(F, F, F, T, T, F, F, F, F, F, T, T, F, F),
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  ) |>
  # hierarquia
  column_spec(
    6,
    underline = c(F, F, F, T, T, F, F, F, F, F, T, T, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Acurácia dos métodos de ML de reconciliação, estratégia *rolling forecast* + validação cruzada com origem móvel"
#| label: tbl-tourism-monthly-results-ml-rolling-fcv

rbind(
  readRDS("data/tourism_monthly/preds_ml/preds/rolling_forecast/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/tourism_monthly/preds_ml/preds/rolling_forecast/resumo_fcv.RDS")[["mase"]]
 ) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped")) |>
  # zone
  column_spec(
    4,
    underline = c(F, F, F, F, F, F, F, F, F, F, T, F, F, F)
  ) |>
  # region
  column_spec(
    5,
    underline = c(F, F, F, T, F, F, F, F, F, F, T, F, F, F)
  ) |>
  # hierarquia
  column_spec(
    6,
    underline = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Acurácia dos métodos de ML de reconciliação, estratégia *fitted values* + validação cruzada aleatória"
#| label: tbl-tourism-monthly-results-ml-fitted

rbind(
  readRDS("data/tourism_monthly/preds_ml/preds/fitted_base/resumo.RDS")[["rmsse"]],
  readRDS("data/tourism_monthly/preds_ml/preds/fitted_base/resumo.RDS")[["mase"]]
 ) |>
 # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped")) |>
  # state
  column_spec(
    3,
    underline = c(F, F, F, F, F, F, F, F, F, F, F, T, F, F),
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  ) |>
  # zone
  column_spec(
    4,
    underline = c(F, F, F, T, F, F, F, F, T, F, T, F, F, F),
    bold = c(F, F, F, F, T, F, F, F, F, F, F, T, F, F)
  ) |>
  # region
  column_spec(
    5,
    bold = c(F, F, F, F, T, F, F, F, F, F, F, T, F, F),
    underline = c(F, T, F, T, F, F, F, F, T, F, T, F, F, F)
  ) |>
  # hierarquia
  column_spec(
    6,
    underline = c(F, F, F, T, F, F, F, F, T, F, T, F, F, F),
    bold = c(F, F, F, F, T, F, F, F, F, F, F, T, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Acurácia dos métodos de ML de reconciliação, estratégia *fitted values* + validação cruzada com origem móvel"
#| label: tbl-tourism-monthly-results-ml-fitted-fcv

rbind(
  readRDS("data/tourism_monthly/preds_ml/preds/fitted_base/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/tourism_monthly/preds_ml/preds/fitted_base/resumo_fcv.RDS")[["mase"]]
 ) |>
 # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped")) |>
  # zone
  column_spec(
    4,
    underline = c(F, F, F, T, F, F, F, F, F, F, T, F, F, F),
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F)
  ) |>
  # region
  column_spec(
    5,
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F),
    underline = c(F, F, F, T, T, F, F, F, F, F, T, T, F, F)
  ) |>
  # hierarquia
  column_spec(
    6,
    bold = c(F, F, F, F, F, F, F, F, F, F, F, F, F, F),
    underline = c(F, F, F, T, F, F, F, F, F, F, T, F, F, F)
  )
```

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Acurácia dos métodos de ML de reconciliação, estratégia *reduced fitted values* + validação cruzada aleatória"
#| label: tbl-tourism-monthly-results-ml-refit

rbind(
  readRDS("data/tourism_monthly/preds_ml/preds/one-step-ahead/resumo.RDS")[["rmsse"]],
  readRDS("data/tourism_monthly/preds_ml/preds/one-step-ahead/resumo.RDS")[["mase"]]
 ) |>
  # aplicando notação científica quando necessário
  sapply(function(x) prettyNum(x, digits = 4)) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped"))
```

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Acurácia dos métodos de ML de reconciliação, estratégia *reduced fitted values* + validação cruzada com origem móvel"
#| label: tbl-tourism-monthly-results-ml-refit-fcv

rbind(
  readRDS("data/tourism_monthly/preds_ml/preds/one-step-ahead/resumo_fcv.RDS")[["rmsse"]],
  readRDS("data/tourism_monthly/preds_ml/preds/one-step-ahead/resumo_fcv.RDS")[["mase"]]
 ) |>
  # aplicando notação científica quando necessário
  sapply(function(x) prettyNum(x, digits = 4)) |>
  # substituindo "glmnet" por "elastic net"
  transform(
    modelo = ifelse(modelo == "glmnet", "elastic net", modelo)
  ) |>
  # substituindo "ranger" por "random forest"
  transform(
    modelo = ifelse(modelo == "ranger", "random forest", modelo)
  ) |>
  kbl(booktabs = TRUE, format = "latex", digits = 3) |>
  pack_rows("RMSSE", 1, 7) |>
  pack_rows("MASE", 8, 14) |>
  kable_styling(latex_options = c("striped"))
```

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Tempo de processamento dos métodos de ML (em horas). Validação cruzada aleatória"
#| label: tbl-tourism-monthly-tempo-ml

# obtendo tempo de processamento
tempo = lapply(c("fitted_base", "rolling_forecast", "one-step-ahead"), function(tipo) {
  lapply(c("xgb", "ranger", "glmnet", "lasso", "ridge", "svm", "lightgbm"), function(learner) {
    preds = readRDS(paste0("data/tourism_monthly/preds_ml/preds/", tipo, "/preds_", learner, ".RDS"))[[2]]
    return(preds)
  })
})

# juntando predições em um único dataframe
tempo = do.call(rbind, tempo)

# convertendo colunas para numérico
tempo = sapply(as.data.frame(tempo), function(x) as.numeric(x))

# nomeando df
colnames(tempo) = c("xgb", "random forest", "elastic net", "lasso", "ridge", "svm", "lightgbm")
rownames(tempo) = c("fitted values", "rolling forecast", "reduced fitted values")

# tabela
tempo |>
  kbl(booktabs = TRUE, digits = 3) |>
  kable_styling(latex_options = c("striped"))
```

```{r}
#| tbl-cap: "Resultados para *visitor nights*: Tempo de processamento dos métodos de ML (em horas). Validação cruzada com origem móvel"
#| label: tbl-tourism-monthly-tempo-ml-fcv

# obtendo tempo de processamento
tempo = lapply(c("fitted_base", "rolling_forecast", "one-step-ahead"), function(tipo) {
  lapply(c("xgb", "ranger", "glmnet", "lasso", "ridge", "svm", "lightgbm"), function(learner) {
    preds = readRDS(paste0("data/tourism_monthly/preds_ml/preds/", tipo, "/preds_", learner, "_fcv.RDS"))[[2]]
    return(preds)
  })
})

# juntando predições em um único dataframe
tempo = do.call(rbind, tempo)

# convertendo colunas para numérico
tempo = sapply(as.data.frame(tempo), function(x) as.numeric(x))

# nomeando df
colnames(tempo) = c("xgb", "random forest", "elastic net", "lasso", "ridge", "svm", "lightgbm")
rownames(tempo) = c("fitted values", "rolling forecast", "reduced fitted values")

# tabela
tempo |>
  kbl(booktabs = TRUE, digits = 3) |>
  kable_styling(latex_options = c("striped"))
```

<!-- ! BUG: Parece ter um erro nos resultados de Spiliotis. Note que tanto em meus resultados de ML (implementados manualmente) quanto os analíticos (utilizando a implementação do {fabletools}, RMSSE < MASE em ambos Tourism (que conferem com FPP3) e Tourism Monthly. Mas, em Spiliotis, RMSSE > MASE tanto para ML quanto os analíticos, indicando que possa haver problemas na execução ou nos dados. -->

# CONCLUSÕES

## Conclusões da dissertação

O rápido crescimento de métodos computacionais e de ML tem permitido a utilização de modelos mais complexos e flexíveis para a previsão hierárquica. Nesse sentido, essa dissertação buscou investigar a aplicação de diversos métodos de ML para a tarefa de reconciliação de previsões incoerentes, com o objetivo de melhorar a acurácia das previsões em diferentes níveis de agregação. Em adição à exploração desses diferentes métodos, a dissertação também contribuiu com elementos para a metodologia de combinação ótima, oferecendo estratégias para aquisição do conjunto de treinamento, bem como a utilização  de validação cruzada com origem móvel na otimização dos hiperparâmetros dos modelos de ML. Essas contribuições foram colocadas a teste em três estudos de caso, um em economia bancária e os demais em bases de dados amplamente utilizadas na literatura de referência, onde puderam ser observado o sucesso das metodologias propostas.

Seguindo um capítulo introdutório, o capítulo 2 revisou os fundamentos de séries temporais hierárquicas e agrupadas, detalhando seus conceitos e álgebra matricial. Foram cobertas as principais abordagens de previsão hierárquica de nível único: a *bottom-up*, as três principais técnicas de *top-down* e o conceito de *middle-out*. Em seguida, foram apresentados os conceitos de coerência e reconciliação a partir da lógica algébrica, assim como a dedução do método de reconciliação ótima por meio de minimização do traço.

As principais contribuições da dissertação se encontram no capítulo 3. Elas se conectam com a metodologia de reconciliação por combinação ótima discutidas na literatura de referência, @spiliotis_hierarchical_2021. No capítulo, são propostas duas estratégias inéditas para aquisição do conjunto de treinamento dos modelos de ML, denominados de *fitted values* e *reduced fitted values*. Essas estratégias podem favorecer a metodologia tanto em relação à acurácia das previsões reconciliadas quanto ao tempo de processamento, a depender da amostra. Além disso, dadas determinadas circunstâncias, argumenta-se a favor da utilização de validação cruzada com origem móvel na otimização dos hiperparâmetros dos modelos de ML.

No capítulo 4, é realizado um estudo de caso no tópico de economia bancária, com dados de empréstimos e financiamentos do Banco do Estado do Espírito Santo. No estudo, foi possível perceber ganhos de acurácia com a metodologia proposta no capítulo 3, com a estratégia *reduced fitted values* obtendo a melhor performance no nível agregado e sua combinação com a validação cruzada com origem móvel gerando a melhor performance nos níveis de meso e microrregião.

Nesse estudo de caso, foram empregados, além do XGB e RF, já utilizados na literatura de combinação ótima, os métodos de ML *elastic net*, em suas configurações *lasso*, *ridge* e *elastic net*, e os métodos SVM e LGBM. Os resultados demonstraram que a aplicação de métodos de reconciliação ótima aumentaram a precisão das previsões em quase todos os níveis hierárquicos, com exceção do nível municipal. Esse resultado vai ao encontro da literatura e já havia sido apontado em @athanasopoulos_forecast_2023, que chama atenção para o fato de que na maior parte dos trabalhos empíricos não é observada melhoria em *todos* os níveis hierárquicos. Não obstante, os ganhos nos demais níveis hierárquicos foram significativos, especialmente nos níveis mais agregados, com destaque para os métodos de ML propostos neste trabalho. Ao contrário dos resultados de @spiliotis_hierarchical_2021, os métodos baseados em árvore de decisão (XGB, RF e LGBM) não superaram o MinT em nenhum nível hierárquico. Já os métodos de regressão regularizada, obtiveram resultados expressivos, com o *elastic net* na estratégia *reduced fitted values* + validação cruzada aleatória alcançando 91% de ganho de acurácia no nível agregado em relação ao método analítico de melhor performance, enquanto o *Lasso* na estratégia *reduced fitted values* + validação cruzada com origem móvel foi o método de melhor performance nos níveis de mesorregião e microrregião, também superando o MinT no nível agregado. Nos demais níveis (agência, verbete e *bottom*), o MinT foi o método de melhor performance.

No capítulo 5, com o objetivo de permitir a comparação com trabalhos anteriores, foram realizados dois estudos de caso adicionais, um com a base de dados *overnight trips* e outro com a base de dados *visitor nights*. Para a primeira, o método SVM na abordagem *rolling forecast* + validação cruzada aleatória obteve a melhor performance em 3 dos 5 níveis, com os dois níveis restantes liderados pelo *ridge* na estratégia *fitted values* + validação cruzada aleatória e MinT. Já na base *visitor nights*, a estratégia *fitted values* + validação cruzada aleatória foi a que obteve os melhores resultados, quando aplicada ao método *elastic net* em sua configuração *ridge*.

## Pesquisa futura

Coincidentemente, cada estudo de caso dentre os três apresentados contaram com uma combinação de método e estratégia diferente como metodologia vencedora. Para pesquisas futuras, pode-se investigar se a performance dos diferentes métodos e estratégias estão relacionadas às características das séries temporais, por exemplo:

- Existem características específicas das séries temporais que favorecem o uso de métodos de ML em níveis agregados ou desagregados?
- Existe motivo teórico para que os métodos de ML baseados em árvore de decisão não tenham superado os demais em nenhum nível hierárquico?
- Os efeitos do ruído de previsão: se os diferentes métodos e estratégias exibem aumento ou deterioração de performance quando as previsões individuais são mais ou menos ruidosas (i.e. se a variância do erro das previsões individuais é maior ou menor).
- Os efeitos de correlação entre as séries: se os métodos e estratégias exibem aumento ou deterioração de performance quando as séries temporais no menor nível hierárquico são mais ou menos correlacionadas.
- Os efeitos de componentes sazonais: se os métodos e estratégias exibem aumento ou deterioração de performance quando as séries temporais do menor nível hierárquico possuem ou não componentes sazonais.
- Os efeitos do tamanho da hierarquia: verificar se os métodos e estratégias exibem aumento ou deterioração de performance quando a hierarquia é mais ou menos profunda (i.e. possui mais ou menos níveis hierárquicos).