Hierarchical Time Series Forecasting of Fire Spots in Brazil: A Comprehensive Approach

: This study compares reconciliation techniques and base forecast methods to forecast a hierarchical time series of the number of fire spots in Brazil between 2011 and 2022. A three-level hierarchical time series was considered, comprising fire spots in Brazil, disaggregated by biome, and further disaggregated by the municipality. The autoregressive integrated moving average (ARIMA), the exponential smoothing (ETS), and the Prophet models were tested for baseline forecasts, and nine reconciliation approaches, including top-down, bottom-up, middle-out, and optimal combination methods, were considered to ensure coherence in the forecasts. Due to the need for transformation to ensure positive forecasts, two data transformations were considered: the logarithm of the number of fire spots plus one and the square root of the number of fire spots plus 0.5. To assess forecast accuracy, the data were split into training data for estimating model parameters and test data for evaluating forecast accuracy. The results show that the ARIMA model with the logarithmic transformation provides overall better forecast accuracy. The BU, MinT(s), and WLS(v) yielded the best results among the reconciliation techniques.


Introduction
Forest fires are a recurring problem in Brazil, with 222,798 fire spots recorded nationwide in 2020 [1].These wildfires significantly impact the ecosystem, infrastructure, and population.The Brazilian biomes host a rich biodiversity that can be partially lost due to these fires.Amazon forest fires substantially affect air quality and public health [2].The quality and quantity of water can also be affected if the fires reach the vegetation around watersheds.Fires can impact crops and, in the long term, compromise food production if soil fertility is affected by the flames [1].Producing statistics regarding fire spots is crucial for implementing more efficient control measures.Forecasting forest fires is pivotal in designing effective policies for preventing and controlling these phenomena [3].
The number of fire spots can be expressed as a time series, whether at the municipal, biome, state, or national level.Several studies have analyzed and forecast time series related to the number of fire spots in Brazil.For example, ref. [4] used multiple linear regression and non-seasonal ARIMA techniques associated with meteorological variables to model and forecast forest fires in the Pantanal of Mato Grosso do Sul in Corumbá.The logarithm of the number of fire spots was introduced as the dependent variable, and maximum temperature, relative humidity, and solar radiation were introduced as independent variables.The results showed that the ARIMA model was more suitable for modeling and predicting forest fires [4].Another study, conducted in the same location, used artificial neural networks (ANNs) to forecast the number of fire spots and showed that it was possible to predict the number of fire spots with a predictive power of 84.8% using a set of meteorological variables as predictors, and 99.4% using only the time series of the number of fire spots as a predictor [5].However, existing studies on fire spot forecasting have been conducted for specific biomes and/or regions, with no study considering making forecasts for all municipalities and biomes in Brazil.
Moreover, the number of fire spots per municipality, biome, and Brazilian territory can be organized hierarchically to form a hierarchical time series (HTS), which refers to a time series that can be aggregated at different levels based on products, geography, or other characteristics [6][7][8].An example of a hierarchical time series is the sales of a particular item in a country, where total sales can be disaggregated by state, and each state's sales can further be disaggregated by municipality.Forecasting hierarchical time series is typically done for all levels, and they should be "coherent", meaning that aggregate forecasts should be equal to the sum of the corresponding disaggregated forecasts [9].For instance, the sales forecast for a product in municipalities should result in state-level sales, and the sum of state-level forecasts should equal the total sales forecast for the country.
The first methods for forecasting reconciliation in hierarchical time series (HTS) were the top-down, the bottom-up, and the middle-out.These involve generating forecasts for one aggregation level and then using them to create coherent forecasts for the other levels.The bottom-up approach requires forecasts for all series at the lowest level, which are then aggregated to obtain higher-level forecasts.The top-down method requires a forecast only for the total series, which will be disaggregated to obtain forecasts for the other series at lower levels.The middle-out approach combines the bottom-up and top-down approaches to make forecasts.There is no consensus in the literature on which method is more efficient, as it may depend on the number of series in the hierarchy and the data quality [10].
More recently, methods have emerged that use all aggregation levels to generate coherent forecasts.These methods combine forecasts for all levels to provide the best aggregate forecasts [6] and are known as optimal combination or optimal reconciliation [10,11].This approach requires estimating the variance and covariance matrix of forecast errors, leading to five variations of this method [9].
To generate coherent forecasts, individual forecasts of the series in the hierarchy are needed without considering any constraints; these are called incoherent forecasts or base forecasts.This paper considers three widely used methods for generating the base forecasts.Exponential smoothing, known as ETS [12][13][14], forecasts the time series as a weighted average of past observations, with weights exponentially decreasing as observations move away from the end of the series.The autoregressive integrated moving average (ARIMA) model generates forecasts based on the autocorrelation of the data.The Prophet, created by Facebook, generates forecasts considering a trend component, seasonality, a component capturing holiday effects, and an error term [15].
Therefore, this study aims to forecast the number of fire spots for each Brazilian municipality, biome, and whole territory by considering different reconciliation techniques and base forecast methods in a three-level hierarchical time series structure.The number of fire spots between January 2011 and December 2022 was obtained from the Brazilian National Institute for Space Research (INPE) and organized by [16].The autoregressive integrated moving average (ARIMA), the exponential smoothing (ETS), and the Prophet models were tested for baseline forecasts, and nine reconciliation approaches, including top-down, bottom-up, middle-out, and optimal combination methods, were considered to ensure coherence in the forecasts.
The rest of this paper is organized as follows.Section 2 describes the data, the methodology behind hierarchical time series and reconciliation techniques, the models used to obtain the base forecasts, and the accuracy measures used to compare the competing models.Section 3 presents a descriptive and exploratory analysis of the data, followed by the results of the accuracy measures for each model.Finally, Section 4 presents the concluding remarks of the study.

The Data
The data used in this study were provided by the QUEIMADAS program of the Brazilian National Institute for Space Research (INPE).The dataset contains the number of fire spots per month in each municipality in Brazil from January 2011 to December 2022 and was organized and made available by [16].The Brazilian territory is divided into six biomes: (i) Amazônia, encompassing approximately 60% of the world's largest rainforest, rich in mineral reserves, and providing 20% of the world's water supply; (ii) Caatinga, featuring a semi-arid climate with remarkable biological diversity and unique species; (iii) Cerrado, acknowledged as the world's richest savanna in terms of biodiversity, which remained largely unchanged until the 1950s when the federal capital was moved to Brasília; (iv) Mata Atlântica, situated along the Brazilian coast and considered the most endangered biome in the country, with only 27% of the original forest cover still preserved; (v) Pampas, characterized by a rainy climate without a dry period and negative temperatures during winter; and (vi) Pantanal, recognized as the planet's most extensive continuous floodplain.
The Brazilian municipalities composed of more than one biome were associated with the predominant biome based on the results of a Brazilian Agricultural Research Corporation (Embrapa) project [17].Figure 1 shows the locations of each of the six biomes in the Brazilian territory.

Hierarchical Time Series
A hierarchical time series is a time series that can be disaggregated into various attributes and different levels.An example of a hierarchical time series is the sales of a specific item in a country, where total sales can be disaggregated by state, and municipalities can further disaggregate each state's sales.Figure 2 illustrates a hierarchical time series with three levels: Total, which is disaggregated into two time series, A and B, and these are further disaggregated into two (AA and AB) and three (BA, BB, and BC) time series, respectively.In a hierarchical time series, obtaining observations at time t is possible by summing the observations from the level below.Considering the series in Figure 2, we have the following equations, , that simplifies as: where y is an n-dimensional vector of all observations at time t, S is the n × m sum matrix, and b t is an m-dimensional vector of all observations at the lower level of the series at time t.

Forecasting Approaches for Hierarchical Time Series
It is possible to represent all forecasting methods using a single matrix notation: where y h represents the coherent forecasts, S is the sum matrix, y h is the base forecasts for all time series in the hierarchy, and G is the matrix mapping the base forecasts, varying in each approach.The following sections briefly describe the different approaches to obtaining coherent forecasts for a hierarchical time series that will be considered in this paper.

Bottom-Up (BU) Approach
The bottom-up (BU) approach involves making forecasts for each time series at the lowest level and then aggregating/summing them to obtain forecasts for the higher levels.The advantage of this method lies in utilizing all available information because the forecasts are made at the lowest level.However, data at the lowest level can be quite noisy, leading to inaccurate forecasts.We can represent this approach using the general form of Equation ( 4), as follows: where 0 m×(n−m) is a null matrix and I m is the identity matrix.The matrix G, mapping the base forecasts, is an identity matrix, as the individual forecasts are directly used as coherent forecasts for the higher levels.The bottom-up approach is useful when the lower levels of the hierarchy significantly contribute to the higher levels or when the data quality at the lower levels is more reliable.

Top-Down (TD) Approach
The top-down (TD) method, or top-down forecasting, involves forecasting the total series and then disaggregating it to the lower levels.For this purpose, proportions p 1 , p 2 , . . ., p m are used to determine how these forecasts will be distributed to obtain the forecast for the lower series.This approach has the advantage of simplicity, as it uses only one series to generate forecasts, which also results in significant information loss [10].In this approach, the matrix G in Equation (3) takes the following form: where p = [p 1 , p 2 , . . ., p m ] is a set of proportions.There are different approaches to obtaining these proportions; ref. [18] propose two: average historical proportions and proportions of the historical averages.The top-down (TD) approach in forecasting hierarchical time series is useful when there is insufficient or unreliable data at the lower levels of the hierarchy, making it challenging to model individual series accurately or when straightforward and computationally efficient forecasting methods are preferred.
The average historical proportions are computed as follows.
where p j represents the historical average proportions of the lower series y j,t in relation to y t .In the results of this paper, this approach is denoted as TDAP.
The proportions of the historical averages are obtained as follows.
where p j represents the historical average value of the lower-level series y j,t in relation to the mean of y t .In the results of this paper, this approach is denoted as TDPA.
The methods described by Gross and Sohl take into account historical proportions but do not consider how these proportions may change over time, which can impact the accuracy of forecasts for lower levels.To address this, ref. [19] proposed the "Forecast Proportions", where forecasts for the lower levels are used to calculate the proportions to be applied.The proportions are calculated as follows: where yj, h (ℓ) is the h-step ahead forecast of the series corresponding to the node that is ℓ levels above j, and Sj, h (ℓ+1) is the sum of the base forecast h steps ahead below the node that is ℓ levels above j, which is directly connected to that node.This approach was denoted as TDFP.

Middle-Out (MO) Approach
In the middle-out approach, or intermediate approach, an intermediate level is chosen, and forecasts are generated for all series at that level.Using these forecasts, the bottom-up approach is applied to the levels above the intermediate level, and the top-down approach is applied to the levels below the intermediate level.This strategy aims to strike a balance between the detailed insights provided by the bottom-up method and the high-level overview offered by the top-down method.The choice of the intermediate level is crucial and may vary depending on the specific characteristics of the time series data and the hierarchy structure.

Minimum Trace (MinT) Optimal Reconciliation Approach
The minimum trace (MinT) optimal reconciliation approach optimally combines base forecasts to produce coherent forecasts.Unlike other approaches, it utilizes all information in the hierarchy to generate unbiased forecasts at all levels of the hierarchy [19].It involves finding the matrix G that minimizes the variances of errors in coherent forecasts y h , as demonstrated by [9].The matrix of variance and covariance of errors in coherent forecasts h steps ahead is given by: where Wh = Var[yT + h − y h ] is the variance-covariance matrix of base forecast errors.Minimizing the variances of errors in coherent forecasts is equivalent to minimizing the trace of the matrix V h , hence the name minimum trace.The matrix G that minimizes the trace of the matrix To obtain y h , it is necessary to estimate the matrix W h .Several approximations have been proposed, giving rise to various approaches.The five main ones are listed below: 1.
W h = k h I for k h > 0, transforms MinT into an ordinary least squares (OLS) estimator and assumes that the matrix G is independent of the data, facilitating calculations [6].However, it assumes independence, leading to information loss.

2.
W h = k h diag( W 1 ) for all h, where k h > 0 and where e t is an n-dimensional vector of residuals from the model generating base forecasts.This approach scales base forecasts using the variance of residuals.In this case, W h can be described as the weighted least squares (WLS) estimator [20].In the results, this approach is denoted as WLS(v).

3.
W h = k h Λ for all h, where k h > 0, Λ = diag(S1), and 1 is a unit vector of dimension m.This implies that lower-level base forecast errors have variance k h and are uncorrelated between nodes.This estimator depends only on the aggregation structure and not on the data, hence termed structural scaling [21].It is denoted by WLS(s) in the results.4.
W h = k h W 1 for all h, where k h > 0. In this approach, error covariance matrices are considered proportional to each other, and the covariance matrix of a one-step-ahead W 1 is estimated using sample covariance.Despite being simple to obtain, it may not be a good estimate when m > t [9].This approach is denoted by MinT(Cov).5.
with W 1,D being a diagonal matrix composed of the diagonal of W 1 , and λ D is the shrinkage intensity parameter.This estimator aims to reduce the sample covariance to a diagonal matrix.λ D is estimated through sample correlation [22], given by: where r ij is the ij-th element of R 1 , the sample correlation matrix one step ahead.This approach is denoted by MinT(s).
The minimum trace (MinT) optimal reconciliation approach aims to minimize the discrepancy between aggregated and coherent forecasts at different hierarchical levels.It achieves this by determining adjustment proportions that minimize the trace of the residual covariance matrix.This approach seeks to optimize the consistency of forecasts throughout the hierarchy, enhancing their accuracy and alignment.In practice, MinT is used as a reconciliation technique to adjust forecasts obtained from different methods, ensuring coherence and adherence to the hierarchical structure.This technique aims to improve the quality of forecasts, especially when distinct methods are applied at different levels of the time series hierarchy.The MinT approach uses all the available information in the hierarchy, but its disadvantage lies in the difficulty of estimating the covariance matrix W h , as it must be positive definite to ensure that its inverse exists, which is not always the case.

Univariate Time Series Forecasting
Individual or base forecasts of the series are required to generate coherent forecasts.This study will consider three models: autoregressive integrated moving average (ARIMA), the ETS class of exponential smoothing models, and Prophet.A brief explanation of each model follows.
The ARIMA model is a combination of three components: autoregressive (AR), moving average (MA), and integrated (I).The autoregressive component involves regressing the variable of interest using its past values.The moving average component for a stationary time series describes the variable of interest as a linear combination of its past forecast errors.The integrated part refers to the difference between consecutive observations of a variable, aiming to make a time series stationary [23].The seasonal ARIMA model can be written as follows: and is denoted by ARI MA(p, d, q)(P, D, Q) s , where • B is the backward shift operator, defined by Bz t = z t−1 , e.g., is the mean average operator of order q; Ps is the seasonal autoregressive operator of order P; is the referencing operator, with d the number of differences to make the time series stationary; • ∇ D s = (1 − B s ) D is the seasonal differences operator; • a t is the white noise.
The selection of the ARIMA model can be done automatically in the R software, version 4.4.0, using the ARIMA function from the Fable package, which uses a variation of the Hyndman-Khandakar algorithm [24].For model selection, the unit root test, minimization of the corrected Akaike information criterion (AICc), bias-corrected version of the AIC for a small sample [25], and maximum likelihood estimation (MLE) are used.
The ARIMA model is a combination of three components: autoregressive (AR), moving average (MA), and integrated (I).The autoregressive component involves regressing the variable of interest using its past values.The moving average component for a stationary time series describes the variable of interest as a linear combination of its past forecast errors.The integrated part refers to the difference between consecutive observations of a variable, aiming to make a time series stationary [23].The ARIMA model can be written as follows: where y ′ t is the differenced time series.This model is called the ARIMA(p,d,q) model, where p is the order of the auto-regressive part, d is the degree of differencing, and q is the order of the moving average part.For data with seasonality, seasonal terms are included in Equation ( 6), and the model would be of the form ARI MA(p, d, q)(P, D, Q) m , where m is the seasonal period.Uppercase letters represent the model's seasonal part, and lowercase letters represent the non-seasonal part.The seasonal term is similar to the non-seasonal term but involves lags of the seasonal period.In the model, the seasonal terms are multiplied by the non-seasonal terms.
The selection of the ARIMA model can be done automatically in the R software using the ARIMA function from the Fable package, which uses a variation of the Hyndman-Khandakar algorithm [24].For model selection, the unit root test, minimization of the Akaike information criterion (AIC), and maximum likelihood estimation (MLE) are used.
In exponential smoothing models, the time series forecasts are a weighted average of past observations, where the weights decay exponentially as the observations move away from the end of the series.The simplest exponential smoothing method is called simple exponential smoothing for data without trend and seasonality.In this method, the forecasts are a weighted average of the data, and the weights decrease exponentially as follows: where 0 ≤ α ≤ 1 is the smoothing parameter.An extension of simple exponential smoothing that allows forecasting data with a trend was created by [13], called the Holt linear model, which adds a trend parameter.Another version of this model is the damped trend, which also includes a parameter that dampens the trend to a flat line in the future [26].A version of the linear Holt model for data with seasonality was proposed by [12,13], and it is known as the Holt-Winters model.This model has two variations: the additive and multiplicative, depending on the nature of the seasonal component.
These models were put into state space form by [27].Each exponential smoothing model consists of an equation describing the observed data and state equations describing how the level (ℓ t ), trend (b t ), and seasonality (s t ) components change over time.Each method described above has two models, one with additive errors and the other with multiplicative errors.Thus, each state space model is denoted as ETS(•, •, •) for (error, trend, seasonality).The selection of the ETS model is done automatically in the R software using the ETS function from the Fable package.
The Prophet model is a recent model created by Facebook [15].The model is ideal for data with strong seasonality and large series.It is based on an additive model composed of four components, which are combined as follows: where g(t) is a trend function, s(t) describes various seasonal patterns, h(t) captures the effects of holidays, and ϵ t is the error term.The model can be selected automatically in the R software using the prophet function from the fable.prophetpackage.

Accuracy Measures
To assess the accuracy of the forecasts, the data were divided into two sets: training data to estimate the model parameters and test data to evaluate its accuracy.The years 2011 to 2021 were used for training, and 2022 was used for testing.This way, 12 forecasts will be obtained for each series to test the model.Two metrics using forecasting errors were employed to measure accuracy.The first is the root mean squared error (RMSE).This scale-dependent measure cannot be used on series with different units.RMSE is defined as follows: where h is the number of forecasts, y t is the actual series values, and y t is the corresponding forecast.The mean absolute scaled error (MASE) and the root mean squared scaled error (RMSSE) were also used, which are useful for comparing series with different units, as they do not depend on the scale of the data [28].The MASE and RMSSE are given by: where n is the length of the time series and m the seasonal period.MASE < 1 and RMSSE < 1 indicate that the proposed method, on average, presents smaller errors than those of a single step of the seasonal naive method.
It is also important to analyze whether the forecasts using different reconciliation approaches yield better results than the base forecasts.The relative root mean squared error (AvgRelRMSE) recommended by [29] can be defined through the geometric mean of the ratios of the mean absolute errors.The AvgRelRMSE is calculated as follows: where RMSE Base i is the RMSE of the base forecasts for series i, and RMSE Rec i is the RMSE obtained for the same series i after reconciliation, and m is the number of time series.One advantage of this metric is that the calculation of (1 − AveRelRMSE) × 100% represents the percentage improvement of each reconciliation approach compared to the base forecasts.In an analogous manner to the AveRelRMSE, the AveRelMASE and AveRelRMSSE can also be calculated.

Data Transformation
To ensure that the base forecasts are positive, data transformation is necessary.This work studied two transformations: log(y + 1) and a square root transformation applied to y + 0.5, where y represents the time series data.In this case, it is necessary to reverse the transformation to obtain the forecasts in the original scale.
In the database used in this study, 140 municipalities had no fire spots, 96 municipalities had only one fire spot, and 13 municipalities had the same distribution of fire spots as another municipality.This makes it impossible to calculate the inverse of the W h matrix of variance and covariance of the base forecast errors for the MinT(s), MinT(c), and WLS(s) approaches, as seen in Section 2.3.4.To solve this problem, random Gaussian noise with zero mean and a standard deviation of 0.001 was added to these municipalities' fire spot time series.After the inclusion of the noise, it was still not possible to calculate the forecasts using the MinT(c) approach.

Flowchart of Methodology
To better understand the methodology used in this study, a flowchart is presented in Figure 3.

Fire Spots in Brazilian Biomes and Municipalities
A fire spot, also known as a heat spot, is a point with a temperature above 47 • C that may indicate fires or burns.Moreover, flames can originate from one or more locations.A spot is detected by monitoring satellites at approximately 700 to 900 km altitude.
In this study, a hierarchical time series with three levels was considered.Level 0 represents Brazil's total number of fire spots from January 2011 to December 2022.Level 1 consists of fire spots disaggregated by Brazilian biome: Amazon, Caatinga, Cerrado, Atlantic Forest, Pampa, and Pantanal.The lowest level, or level 2, represents biomes disaggregated by municipality in Brazil, totaling 5570 municipalities.Thus, a total of 5577 time series are included in the hierarchy.
Table 1 presents the descriptive measures of the number of fire spots in Brazil, its biomes, and municipalities.The data show high variability in monthly fire spots across all hierarchy levels, particularly in the municipalities.The Amazon and Cerrado biomes exhibit the highest numbers of fire spots, with significant variability between months, while the Pampa biome has the lowest number of fire spots and the least variability in the data.When grouped by their predominant biome, municipalities generally show a prevalence of months with few or no fire spots; 75% of the total observations from municipalities where the Caatinga, Atlantic Forest, and Pampa are the predominant biomes reported zero fire spots.Figure 4 displays the time series of fire spots in Brazil and its biomes, representing the series at levels 0 and 1 of the hierarchy.In 2017, Brazil experienced the month with the highest number of fire spots between 2011 and 2022, with over 60 thousand fire spots.The Amazon had the highest number of fire spots among the six biomes, while the Pampa had the lowest.The sudden increase in the number of spots in the Pantanal in 2020 is noteworthy, a year in which the biome lost approximately 30% of its vegetation [1].The seasonal component in the series is visible in Figure 4, being less explicit for the Pantanal due to the peak in spots in 2020 and the Pampa, which does not exhibit a welldefined pattern.Due to the large number of municipalities (5570), it is unfeasible to visualize the series at the last level of the hierarchy.To better understand the distribution of fire spots and consider the discrepancy in the number of spots between municipalities, analyses were conducted regarding the number of spots per area (km 2 ) from 2011 to 2022 in each municipality (Figures 5-8).Most municipalities showed a low number of spots per area.There is no clear trend of increase or decrease in spots over the years (Figures 5 and 6).In the Pantanal, the peak of fires in 2020 is again visible in this graph, and only two municipalities were responsible for most spots in the biome that year (Figure 6). Figure 6 shows that fire spots are not evenly distributed across the country; some municipalities are responsible for most spots in Brazil.From 2011 to 2022, the country's southern region had few fire spots per km 2 , and the state of Santa Catarina, belonging to this region, had no spots in some years.Municipalities in the states of Maranhão and Tocantins had a high number of spots per km 2 throughout all years.These states have the Cerrado biome as their predominant biome.The analysis of the distribution of fire spots per area in each municipality throughout the months of the year indicates the presence of seasonality in the data (Figures 7 and 8).At the beginning of the year, the number of fire spots per area in municipalities is low or zero, and this number begins to increase in July (Figures 7 and 8).The peak of fire spots occurs between September and October, depending on the municipality's biome.Only in cities where the Pampa is the predominant biome is a pattern less clear to identify (Figure 7).In September, the highest number of fire spots per km 2 occurs throughout the country, and by December, this number decreases considerably (Figure 8).Fires are more frequent in the winter and spring months, from the end of June to the end of December, due to low rainfall levels [30].

Hierarchical Time Series Forecasting
The routines for calculating fire spot forecasts for the twelve months of 2022 were implemented using R software, version 4.3.0.The results in this section pertain to reconciled forecasts using the square root transformation of fire spots plus 0.5 and the logarithmic transformation of fire spots plus 1.Both transformations used the ARIMA model for base forecasts.The ETS model was also used for base forecasts with the square root transformation of fire spots plus 1. Reconciled forecasts using ETS and the logarithmic transformation of fire spots encountered convergence errors, and the adjusted values and forecasts were unrealistic, so they were not presented.
Table 2 shows the root mean squared error (RMSE) for the level 0 series, related to the overall number of fire spots throughout Brazil, and the average RMSEs for the level 1 and 2 series.Values in bold indicate the lowest RMSEs at each hierarchical level.The results of the average accuracy for different reconciliation approaches using the square root plus 0.5 transformation in the data and the ARIMA model for base forecasts are presented in the first part of the table.With this setup, at level 0, the OLS approach exhibited the lowest RMSE.At level 1, the WLS(s) approach yielded better results, followed by the MinT(s) approach.At level 2, the MinT(s) and WLS(s) outperformed the other approaches.When considering the ARIMA model for base forecasts and a logarithmic transformation of the number of fire spots plus 1, at level 0, the MinT(s) approach had the lowest RMSE.At levels 1 and 2, BU, MinT(s), and WLS(v) outperformed the others (Table 2).The logarithmic transformation results outperformed those based on the square root transformation when the ARIMA model obtained the base forecasts, especially for BU, MinT(s), and WLS(v).When considering the square root transformation and the ETS model to obtain the base forecasts, BU showed better results at level 0, and MinT(s) at the other two levels.The average RMSE values were generally higher than those obtained using the ARIMA model to obtain the base forecasts (Table 2).Reconciled forecasts using ETS and the logarithmic transformation encountered convergence errors, and the adjusted values and forecasts were unrealistic, so they were not presented.When considering the Prophet to obtain the base forecasts, the overall best performance was obtained by BU, followed by MinT(s) and WLS(v) for both transformations (Table 2).Table 3 is similar to Table 2, but uses MASE as the accuracy measure.The ARIMA model with logarithmic transformation presented the lowest average MASE values, particularly with the MinT(s) approach for level 0 and WLS(v) for levels 1 and 2. Comparing the average MASEs by hierarchical level, level 0, representing Brazil, showed the best results in terms of accuracy, followed by level 1.The MASE values for municipalities were, on average, greater than 1 for all models and approaches.
Table 4 presents the results for the RMSSE, which are similar to those of RMSE and MASE (Tables 2 and 3) in terms of the best model and approach.Unlike for MASE (Table 3), at level 2, the average RMSSE was less than one for some approaches.The base forecasts without considering any hierarchical structure, especially at levels 0 and 1, resulted in worse outcomes when compared to some reconciliation approaches (Tables 2-4).In addition to calculating the RMSE, MASE, and RMSSE, it is important to know which approaches produce better forecasts than the baseline.If they are not better, hierarchical time series forecasting loses relevance.Tables A1-A3 of Appendix A show the AveRelRMSE, AveRelMASE, and AveRelRMSSE, respectively, for each level of the hierarchy considering the ARIMA, ETS, and Prophet base forecast models and the two transformations.A value of the AvgRelRMSE, AveRelMASE, and AveRelRMSSE equal to one indicates that the reconciled forecast is equal to the base forecast at that level.An example is the bottom-up approach, which uses base forecasts from municipalities, level 2, and sums them to obtain forecasts for the above levels.A value of AveRelRMSE, AveRelMASE, and AveRelRMSSE lower than one indicates that the reconciled forecasts showed improvements compared to the base forecasts, with bold values being the lowest in Tables A1-A3.
Here, we summarize the results for the AveRelRMSE (Table A1 of the Appendix A).The interpretation for the AveRelMASE (Table A2 of the Appendix A) and AveRelRMSSE (Table A3 of the Appendix A) are done similarly.
Considering the ARIMA model and the square root transformation of fire spots plus 0.5, at level 0, only the OLS and MO approaches had a value of less than one, showing an improvement of 11.1% and 5.0%, respectively, in RMSE compared to the base forecasts.At level 1, WLS(s) was the one that showed a higher improvement, with a value of 5.8%.At level 2, TDFP showed an improvement of 5.3%.Considering the ARIMA model and the logarithmic transformation of fire spots plus 1, at level 0, only the forecasts from the BU and WLS(v) approaches did not show improvements in the RMSE compared to the base forecasts.At this level, MinT(s) obtained the best result, with an improvement of 22.2% in the RMSE value.At level 1, out of the nine approaches, four of them showed improvements in the RMSE up to 20.0%.In the last level, no method showed improvement compared to the base forecasts (Table A1).When considering the Prophet model for the base forecasts, the overall best reconciliation method was the BU, but without an RMSE improvement in level two.Overall, the MinT(s) and TDFP approaches outperformed the remaining when considering the logarithm transformation and the ARIMA model to obtain the base forecasts.
The AvgRelRMSE of the approaches using ETS for base forecasts and the square root transformation of fire spots plus 0.5 is also presented in Table A1.At level 0, excluding the top-down approaches that use base forecasts at this level, all approaches showed improvements in RMSE.BU showed an improvement of 76.3% at this level.At level 1, only WLS(s) and TDPA showed improvements in RMSE, and six approaches showed improvements in RMSE at level 2, with TDPA having the highest improvement.It is also noted that the MO and TDPA approaches showed improvements at all levels.
When considering the Prophet model for the base forecasts, the most effective reconciliation method overall was the BU, except for level 0 in the case of the square root transformation, where the TDPA showed superior performance.
Despite the percentage improvements in RMSE, it is important to consider that the reconciled forecasts using each base model and transformation have different magnitudes.For example, the ETS model with the square root transformation had the overall highest RMSE values (Table 2) compared to the prophet with the logarithm transformation.Thus, reconciled forecasts using the ETS model show significant improvements over base forecasts, but they are not better than those using the Prophet model.
Figures 9 and 10 show the boxplots for the RMSE for the nine forecasting reconciliation approaches and the three models and two transformations under consideration for level one (biomes) and level two (municipalities), respectively.For level one (Figure 9), the WLS(s) seems to be the better-performing reconciliation approach, and the ETS model with the square root transformation is the best to obtain the base forecasts.For level two (Figure 10), when not considering the "outliers", a good overall performance was obtained for the TDFP reconciliation approach, and the ETS with the square root transformation seems to be the overall winner in performance.Figures A1 and A2 show the boxplots for the MASE for the nine forecasting reconciliation approaches and the three models and two transformations under consideration for level one (biomes) and level two (municipalities), respectively.Figures A3 and A4 show the boxplots for the RMSSE for the nine forecasting reconciliation approaches and the three models and two transformations under consideration for level one (biomes) and level two (municipalities), respectively.The overall conclusions from Tables A1-A3 are in the same line as those for the RMSE (Figures 9 and 10).

Concluding Remarks
In this study, we comprehensively analyzed various reconciliation approaches to forecasting the number of fire spots, considering different base forecast models and data transformations.The results offer valuable insights into the effectiveness of these approaches at different hierarchical levels, ranging from biomes to municipalities.
Our findings indicate that the choice of the reconciliation method is crucial and depends on factors such as the base forecasting model and the transformation applied to the data.At level 0, the bottom-up (BU) approach demonstrated stronger performance, particularly when using the Prophet model and the logarithm transformation.However, alternative methods, such as WLS, BU, and MinT, exhibited competitive results at other levels.
Moreover, the influence of data transformations, including square root and logarithmic transformations, was evident in the performance of reconciliation approaches.Notably, when combined with the ARIMA or the Prophet model, the logarithmic transformation showcased superior results for certain approaches at various levels.
The AveRelRMSE, AveRelMASE, and AveRelRMSSE metrics provided a nuanced understanding of the improvements achieved by each reconciliation approach compared to baseline forecasts.This information is essential for practical decision-making, emphasizing the importance of not only minimizing forecast errors but also ensuring improvements over the baseline.
It is worth noting that the reconciliation methods exhibited diverse performance across different levels of the hierarchical structure, suggesting the need for adaptive strategies based on the specific forecasting context.
In summary, this study contributes valuable insights into applying forecasting reconciliation in fire spot forecasts.The results presented here guide practitioners and decision-makers in choosing reconciliation methods tailored to their specific forecasting scenarios.As the field continues to evolve, future research may explore additional factors influencing reconciliation effectiveness and expand the applicability of these findings to diverse domains.

Figure 1 .
Figure 1.Geographic localization of each of the six Brazilian biomes.

Figure 3 .
Figure 3. Flowchart of the methodology used in the paper.

Figure 4 .
Figure 4. Time series of the monthly number of forest fire spots in Brazil and each of the six Brazilian biomes.

Figure 5 .
Figure 5. Boxplots for the annual number of fire spots in Brazilian municipalities divided by their area and grouped by the municipality's predominant biome.

Figure 6 .
Figure 6.Maps of the total fire spots in Brazilian municipalities divided by their area from 2011 to 2022.

Figure 7 .
Figure 7. Boxplots for the monthly number of fire spots in Brazilian municipalities divided by their area and grouped by the municipality's predominant biome, considering the years from 2011 to 2022.

Figure 8 .
Figure 8. Maps of the total forest fire spots in Brazilian municipalities divided by their area per month, considering the years from 2011 to 2022.

Figure 9 .
Figure 9. Boxplots for the root mean squared error for 12 ahead forecast of the number of fire spots for all six Brazilian biomes considering the nine forecasting reconciliation approaches, using (i) the square root transformation and the ARIMA model for base forecasts; and (ii) the logarithmic transformation and the ARIMA model for base forecasts; (iii) the square root transformation and the ETS model for base forecasts; (iv) the square root transformation and the Prophet model for base forecasts; and (v) the logarithmic transformation and the Prophet model for base forecasts; red line: mean of the medians.

Figure 10 .
Figure 10.Boxplots for the logarithm of the root mean squared error plus one for 12 ahead forecast of the number of fire spots for all 5570 Brazilian municipalities considering the nine forecasting reconciliation approaches, using (i) the square root transformation and the ARIMA model for base forecasts; and (ii) the logarithmic transformation and the ARIMA model for base forecasts; (iii) the square root transformation and the ETS model for base forecasts; (iv) the square root transformation and the Prophet model for base forecasts; and (v) the logarithmic transformation and the Prophet model for base forecasts red line: mean of the medians.

Figure A3 .
Figure A3.Boxplots for the root mean squared scaled error 12 ahead forecast of the number of fire spots for all six Brazilian biomes considering the nine forecasting reconciliation approaches, using (i) the square root transformation and the ARIMA model for base forecasts; and (ii) the logarithmic transformation and the ARIMA model for base forecasts; (iii) the square root transformation and the ETS model for base forecasts; (iv) the square root transformation and the Prophet model for base forecasts; and (v) the logarithmic transformation and the Prophet model for base forecasts; red line: mean of the medians.

Figure A4 .
Figure A4.Boxplots for the root mean squared scaled error for 12 ahead forecast of the number of fire spots for all 5570 Brazilian municipalities considering the nine forecasting reconciliation approaches, using (i) the square root transformation and the ARIMA model for base forecasts; and (ii) the logarithmic transformation and the ARIMA model for base forecasts; (iii) the square root transformation and the ETS model for base forecasts; (iv) the square root transformation and the Prophet model for base forecasts; and (v) the logarithmic transformation and the Prophet model for base forecasts; red line: mean of the medians.

Table 1 .
Descriptive measures of the number of fire spots in Brazil, its biomes, and municipalities.

Table 2 .
Average RMSE per hierarchical level using (i) the square root transformation and the ARIMA model for base forecasts; (ii) the logarithmic transformation and the ARIMA model for base forecasts; (iii) the square root transformation and the ETS model for base forecasts; (iv) the square root transformation and the Prophet model for base forecasts; and (v) the logarithmic transformation and the Prophet model for base forecasts.

Table 3 .
Average MASE per hierarchical level using (i) the square root transformation and the ARIMA model for base forecasts; (ii) the logarithmic transformation and the ARIMA model for base forecasts; (iii) the square root transformation and the ETS model for base forecasts; (iv) the square root transformation and the Prophet model for base forecasts; and (v) the logarithmic transformation and the Prophet model for base forecasts.

Table 4 .
Average RMSSE per hierarchical level using (i) the square root transformation and the ARIMA model for base forecasts; (ii) the logarithmic transformation and the ARIMA model for base forecasts; (iii) the square root transformation and the ETS model for base forecasts; (iv) the square root transformation and the Prophet model for base forecasts; and (v) the logarithmic transformation and the Prophet model for base forecasts.

Table A3 .
AveRelRMSSE per hierarchical level using (i) the square root transformation and the ARIMA model for base forecasts; (ii) the logarithmic transformation and the ARIMA model for base forecasts; (iii) the square root transformation and the ETS model for base forecasts; (iv) the square root transformation and the Prophet model for base forecasts; and (v) the logarithmic transformation and the Prophet model for base forecasts.