1 Introduction
A variety of economic systems consists of a certain class of subsystems, which form a fixed hierarchical structure, i.e.:

A country considered with respect to monthly Gross National Product with a splitting into regions and subregions

A national balance of trade with a splitting into branches and subbranches of industry, services and agriculture

Weekly total inflows and outflows of current accounts for a certain group of clients (priority, individual, corporate clients, or gender, or demographic groups) considered in 5min consecutive intervals

A turnover of a company with regard to product lines and/or client target groups

Social and health care insurance costs divided into age and place of living segments

Social inequalities of households with regard to education level, religious faith, political choice or ethnicity

Social and health costs associated with environment pollution for a certain region divided into subregions
Both from a theoretical as well as from practical point of view it is especially important to find a reliable method of modeling and forecasting a time evolution of a system similar to the above systems exhibiting a hierarchical structure. The method should be computationally tractable.
The issue is very closely tied with the reconciliation of forecasts  a problem known from the econometric literature (see Shlifer and Wolff (1979), Kohn (1982), Weale (1988), Kahn (1998), and Fliedner (2001)).
It is often observed that forecasts prepared for lower levels of a hierarchy do not sum to forecasts prepared for upper levels and a top level of the system. That fact may be caused for example by an application of different measurement methodologies, different precision for different levels.
In this paper we focus our attention on a problem of forecasting a hierarchic system describing a day and night air pollution in Silesia region in Poland. The region is divided into five subregions.
The day and night air pollution is treated as a realization of a functional random variable. Hence we consider a forecasting of a hierarchical functional time series (HFTS). A predictor possessing good statistical properties in this case is directly connected with an opportunity of designing costeffective proecological regional politics, which optimize a social welfare being a function of a day and night air pollution. It is worth stressing, that while using functional time series (FTS) framework instead of well known onedimensional time series setup, we forecast whole day and night periods instead of predicting hour after hour. Note also, that using FTS one can easily model and predict nonequally spaced time series, which may causes fundamental problems for analysts using classical ARIMA, SARIMA methodology (see Kosiorowski 2016, and Górecki et al. 2016). Our paper critically discusses the best proposals known from the literature (see Shang and Hyndman 2017) and compare them with our proposals. The rest of the paper is organized as follows. Section 2 discusses the concept of hierarchical time series. Section 3 discusses the concept of hierarchical functional time series. Section 4 discusses the methods of HFTS forecasting. Section 5 contains the empirical study – an example of a macromodel for a day and night air pollution in Silesia region.
2 Hierarchical time series
Hierarchical time series (HTS) is a time series, where some fixed, often natural, hierarchy is imposed. In other words, HTS can be considered as a time series, where at each time we have insight to the values for any single variable at any level in the structure with a fixed hierarchy. In Figure 1 an example visualization of hierarchical time series at moment
is depicted. In the Figure 1, the observation made on the top level is divided into two sublevels or levels, and the observation made on the is divided into two levels, but the division might be quite different, and the only constraint is, that any level could be divided into finitely many number of levels and the total number of observations is finite. Obviously, one can compute a forecast for all series at all hierarchy levels independently, but the forecast at the lower level do not sum to the forecast at upper level. Hence, no reconciliation is made. In the hierarchical setup the forecasting might be done in the following manners. The forecast is made on the bottom level of the hierarchy. Subsequently, the aggregation of the obtained forecast, basing on some historical data, is made on the upper level of the hierarchy. This procedure is repeated upwards the hierarchy, until we get a forecast on the top level. The method is called the bottomup method. Conversely we proceed in the topdown method, where the forecast is made on the top level. The disaggregation is made then, so that we obtain a forecast on the lower levels of the hierarchy. The method are often mixed, as we obviously, for some reasons, might be interested in the forecast on some intermediate level of the hierarchy. Then the forecast is aggregated upward the hierarchy, and disaggregated downward the hierarchy. The methods do not take into account the correlation structure of the hierarchy. Prediction intervals for the forecasts are undefined as well. The more detailed discussion and references the interested reader may find in Shlifer and Wolff (1979), Kohn (1982), Weale (1988), Kahn (1998), and Fliedner (2001).In their paper Hyndman et al. (2011) proposed a novel optimal combination forecast method for HTS. Their proposal is based on independently forecasting all series at all levels of the hierarchy and then using a regression model to obtain a reconciliated forecast. The forecast obtained with their method add across the fixed hierarchy. It is also meanunbiased and under some reasonable assumptions has minimum variance among linear combination of independent forecasts. They represent the fixed hierarchical structure in the matrix form. This approach allows for the correlations between the series at each level of the hierarchy. However, they mention that some computational problems may occur. They are connected to the inverse of relatively large, sparse matrices and solution of sparse linear least squares problem. Nevertheless, Hyndman et al. (2011) approach enables to obtain a reconciliated forecast for a considered phenomenon reconciliated with individual forecasts obtained on different levels of hierarchy.
3 Functional hierarchical time series
Functional hierarchical time series is a series which consists of functional data, i.e. we consider a hierarchical dataset of functions instead of real numbers or vectors in
The functional data methods are described in monographies of Ferraty and Vieu (2006), Ramsay et al. (2009), and Horváth and Kokoszka (2012). Some statistical tests have been recently developed for the functional framework as well (e.g. Kosiorowski et al. 2017a).
Notice, that the methods developed for the uni or multivariate HTS and described in Section 3 cannot be directly applied in HFTS setup. Many theoretical problems arises here, but note, that even the order of functions cannot be measured as easily as in the univariate case. The naïve forecast is not convincing as well, because it does not take into account a time dependency. The pointwise average can be easily calculated, but we usually do not know the true distribution on the space, from which our data come from, so we can’t straightly assume that the functional expected value exists, which makes the approach unconvincing as well. For the same reason the pointwise moving average seems out of the question. In their paper Shang and Hyndman (2017) proposed their method of HFTS forecasting. Their approach originates from their previous study (Hyndman et al. 2011), described in Section 2, which takes into account the whole hierarchical structure of the data. The reconciliated forecast for a fixed hierarchical structure with levels takes a following shape
(3.1) 
where denotes an forecast obtained for the functional time series at level and
denotes a certain generalized least squares estimator. The HFTS structure at time
is described by a matrix equation(3.2) 
where vector , that is, it is containing all series at all levels of hierarchy, is a vector representing the series at the lowest level of the fixed hierarchy, and is a finite matrix that shows the connection between the vectors and .
The forecast is made then, that is:
(3.3) 
where is a matrix of forecasts made for all series at all levels of the fixed hierarchy, is an unknown multivariate expected value of a forecast distribution for the most disaggregated series and represents the errors of the reconciliation.
Note that the level forecasts are obtained using nonrobust method, which maps functional time series into one dimensional series of functional component scores (for details see Kosiorowski 2014). Components are estimated in the study by Shang and Hyndman (2017) with the generalized least squares method, i.e.
(3.4) 
where a diagonal matrix estimates variances of series forecasts. The final forecast stems from the equation
(3.5) 
Hyndman and Shang’s method has some important advantages and disadvantages. The Hyndman and Shang forecasts are aggregate consistent – they satisfy an aggregation constrains and are meanunbiased. However, the method is very computationally demanding and sophisticated due to the sparse linear least squares problem, and necessity of computing generalized inverses of large and sparse matrices (see 3.2). The method may be therefore robustified (for details see Kosiorowski et al. 2017b).
Note that Shang and Hyndman’s approach is equipped with an outright internal mechanism of forecasts reconciliation. Our proposal – which is described further – the reconciliation of forecasts is a byproduct of a fact that MBD depth is nontransitive.
Hyndman and Shang reduce the problem of functional data forecasting to functional principal component regression: functions are represented in empirical principal components base, then they use Hyndman’s functional regression basing on onedimensional stationary time series modelling (see Kosiorowski 2014) and the authors assume that residual functions are approximately stationary (we don’t make the restriction).
The last analyzed approach is ours. In a paper Kosiorowski et al. (2017c) we presented our double functional median method and compared our method with Hyndman and Shang method as a reference approach.
Modified band depth (MBD, see LópezPintado and Romo, 2007 and 2009) of curve with respect to functional sample estimates the curves’ frequency of being in the center. Note, that Zuo and Serfling (2000) formulated general conception of statistical depth function and NietoReyes and Battey (2016) have proved that a depth for functional data is correctly defined.
Nevertheless, we have a sample of functions . Firstly we need to define sets of the following form
Consequently, MBD can be defined, a functional depth, which takes into account a proportion of "time", when is in the band made with two functions, i.e.
(3.6) 
Subsequently, the nested regions for the chosen functional depth can be constructed, that is, consider . The median with respect to the considered functional depth is the most central observation. Assume that we have functions, i.e., and
denotes sample functional depth of function with respect to sample . We define a sample median as
(3.7) 
If more than one function is achieving the depth maximum value, the median is defined as the average of the curves maximizing depth. We use then a moving functional median:
(3.8) 
where is a moving window of a length with an end in a moment , that is,
Our method can be described in the following steps:
First step: we calculate the moving functional median related to the MBD or another functional depth for each unit at the lowest level of hierarchy, i.e., In empirical example analyzed in Section 5 for each town and at moment , we compute a functional median from a moving window of length with respect to the functional depth :
Second step: we calculate for the lowest but one level of hierarchy, a functional median from medians calculated in the first step.
We repeat the second step until we calculate the functional median for the top level of the hierarchy.
In our empirical example, the second step is the last one, and finally, we obtain a forecast for
(3.9) 
Hierarchical structure of the data is taken into account in the process of computing functional median of lower level functional medians, as translation of a single functional observation into neighbouring unit alters the outcome (for details see Kosiorowski et al. 2017c).
4 A critical overview of HFTS approaches
In a general case, an uncertainty evaluation of the HFTS forecast is an open issue. Due to insufficient theoretical background for conducting a precise statistical inference, in our approach we decided to expand ideas indicated by LópezPintado et al. (2010).
Shang&Hyndman (2017) has obtained a
representation of functions in the space with a Fourier basis. The Fourier basis is adjusted to the functional data they consider, because the data they analyzed were expected to be periodical.
Afterwards they transformed functional time series into a family of one dimensional principal component scores series.
Then a maximum entropy bootstrap methodology proposed by Vinod and de Lacalle (2009) and implemented in meboot R package, which is appropriate for time series setting, has been used.
Although this simplification of the problem seems to be attractive, it divests an analyst the richness of behaviors a functional time series in comparison to one dimensional time series.
We recommend using functional boxplots and adjusted functional boxplots (one can focus on sizes of boxes and central regions), which realizes an idea of bootstrap for functional time series for the forecast uncertainty rough evaluation (for details see LópezPintado et al. 2010 and Sun and Genton 2011).
It does not make much sense to consider pointwise properties of the considered predictors.
We usually do not know the true distribution on the
space, from which our data come from. Thus even the existence of the functional expected value (mean) cannot be assumed. Hence, we concentrate our attention on the medianunbiasedness. Let’s remind, that an estimate of a onedimensional parameter is medianunbiased if, for fixed parameter value, the median of the distribution of the estimate is at the parameter value, which simply means that the estimate underestimates just as often as it overestimates (Brown, 1947). The classical medianunbiasedness properties have been studied previously (e.g. Pfanzagl, 1970 and 1979). In the functional setting, we choose the proper functional depth and thus we obtain the median induced by the chosen depth. The functional medians induced by popular depth exist for very wide class of processes (in contrary to the functional mean existence). We conclude, that the functional median obtained with respect to the chosen functional depth is intrinsically a medianunbiased estimator. Moreover, our double median method is not only medianunbiased, but also consistent (for details see Gijbels and Nagy 2015, Nagy et al. 2016 and Kosiorowski et al. 2017c).
Shang and Hyndman method (2017) depends on quite effective but nonrobust onedimensional time series methodology applied to series of principal component scores. It depends also on nonrobust dispersion matrix estimator. The matrix is a kind of a design matrix used to obtain a proper forecasts reconciliation. The robustness of our method to outliers does not heavily depend on the type of functional outliers. It is surprising, because we have expected, that it should be different for the functional shape outliers, functional amplitude outliers, and for functional outliers with respect to the covariance structure (e.g. see ArribasGil and Romo 2014 and Tabalerroni 2017). After conducting several simulations (see Kosiorowski et al. 2017c) we have come to the conclusion, that the double median method is more robust. Hyndman and Shang state the opposite in their paper, but note, that they considered a FraimanMuniz depth, while we have considered MBD, that looks like better designed for the considered empirical example.
For a fixed , a volume of the central region may be treated as a dispersion measure (see Liu et al., 1999), and thus comparing functional boxplots is a relevant way to compare "effectiveness" of the considered methods (see figures 3, 4, 9). A comparison of functional time series predictor "effectiveness" may be also conducted in terms of speeds of expansion of central regions treated as a functions of (scale curve, see LópezPintado et al., 2010). This approach is not only nonparametric but also a momentfree dataanalytic method. It imitates the multivariate case and seems to be the best solution in the functional case as well, because assumptions on the datagenerating process are hard to be stated precisely. Remember, that there in no Lebesgue measure analogue in the space.
The double median method of forecasting HFTS, is faster than Hyndman and Shang’s method, moreover, it is less computationally and memory intensive then theirs. Precisely, to compare a computational complexity of both methods, we have considered empirical functional time series related to day and night air pollution monitoring in the selected towns. The monitoring was conducted for 181 days. In other words, in the beginning, we considered dataset consisted of six matrices, each of dimension . For comparing two forecasting methods, we have considered forecasts obtained basing on moving window of length 10 observations. A time of calculation of the forecasts using Shang & Hyndman method was ca 13 min 10 sec, whereas using the proposed double median method was ca 2 min 30 sec (we used DepthProc R package). In both cases, we used the same software and hardware environment (WIN8, Intel Core I7 Mobile, 16 GB RAM). Note that Shang and Hyndman’s (2017) and Hyndman et al. (2011) indicated the inconveniences of their methods, which are related to an application of the generalized least squares applied to big and sparse design matrices. They stated some methods to bypass the inconveniences, but the remedies are insufficient in big data analysis.
5 Empirical study: A day and night PM 10 air pollution in Silesia region
Air pollution consist of different substances, i.a. sulphur dioxide, nitrogen dioxide, ozone, carbon monoxide, benzene, particulate matter PM2,5 and particulate matter PM10  all particles of a diameter 10 micrometers or less. Air pollution has a huge negative impact on people’s health.
Air pollution monitoring is conducted in Silesian Province in Poland.
Measurement is done at a certain number of stations placed in the Region. The organisation responsible for the monitoring is
Wojewódzki Inspektorat Ochrony Środowiska (WIOŚ, Regional Inspectorate of Environmental Protection) in Katowice.
The institution posses 28 measurement stations. We analyze data coming from 5 of that 28 stations in order to present our method, but the number of stations does not limit our method.
Decision maker who has in a disposal measurement from certain number of stations is interested in aggregation of the data.
Note, that the easiest aggregate is an arithmetic mean or a moving arithmetic mean of measurements done in all the stations.
The two aggregates are often used in practice by the local government.
However, simplicity seems to be the only advantage of that method (see Section 3 and our paper Kosiorowski et al. 2017c).
Main goal of the paper is to find the aggregate, to the best fit to the regional policy.
5.1 Empirical dataset under study description
Dataset from WIOŚ website http://powietrze.katowice.wios.gov.pl has been analyzed to illustrate our method.
We have analyzed PM10 concentration in the air for five measurement stations: Gliwice (Gli) with a population of 182,155, Katowice (Kat) with a population of 304,063, DąbrowaGórnicza (Dab) with a population of 121,902, BielskoBiała (Bie) with a population of 172,407 and Częstochowa (Cze) with a population of 227,184 (population data come from 2015).
First three towns are part of "Upper Silesian Urban Area" (its population is about 3 million).
BielskoBiała and Częstochowa are the largest towns of Silesian Region, that are not part of the "Upper Silesian Urban Area". Data comes from the period of 181 days from 1 September 2016 to 28 February 2017. We obtain forecasts for each town, but we are rather interested to obtain a forecast for the whole Silesian Region, and we shall keep in mind that emission of pollution and weather conditions (i.e. landform and windrose) are very different in each town. Moreover, some of that factors, i.e. wind, are time variant, so we should treat the observations as a functions and treat the trajectories as functional data objects.
It is cumbersome in air pollution context to decide, how to compute air pollution on the whole Silesian Region level, as obviously only data from certain stations are available. We decided to compute an aggregate representing air pollution in the Silesian Region to be a weighted average of pollution in each town, where weights are proportional to the town population. This approach is compatible with our assumption, that social cost associated with air pollution, is linearly proportional to town population. The assumption has been applied in double functional median method and in Shang & Hyndman method.
Figure 2 presents the raw data of 181 curves for the analyzed five stations, which show the PM10 concentration in the atmosphere in on vertical axis. Above there is a functional boxplot for all curves ( curves) computed with MBD.
Figure 3 presents PM10 concentration in the air forecast in calculated with moving functional average (moving window equals 10) for five considered stations. Above there is a functional boxplot for the average of all averages for five stations computed every day. The moving average seems to be the easiest rational method of forecasting, which is used by the decision maker.
Figure 4 presents PM10 concentration in the air double functional median forecast in calculated with moving functional median (moving window equals 10) for five stations. Above the forecast for the Silesia region calculated with double moving functional median. The median was calculated with the use of MBD.
Hyndman and Shang (2017), in order to estimate prediction uncertainty, have used maximal entropy bootstrap for time series method proposed by Vinod and de Lacalle (2009), because their method is basing on representation of functional time series as a family of onedimensional time series of functional principal component scores.
In case of our method, in order to estimate prediction uncertainty we have used volumes of central regions (see functional boxplots) implemented in Rpackages fda (Ramsay et al. 2009) and DepthProc (Kosiorowski and Zawadzki, 2017). We have also compared quality of our forecast with the forecast of Shang and Hyndman through the comparison of sum of the differences between the observed curves and of the forecasted curves. We have also compared median absolute deviation (MAD) of the integrated differences between the observed curves and of the forecasted curves.
Table 1 contains a MAD comparison of our forecasts with Shang and Hyndman’s. The functional boxplots can be also used to compare the two methods. Figure 5 presents five functional boxplots for the average sum of the differences between the observed curves and the curves forecasted with the double functional median method.
Above there is a functional boxplot for the forecasts for Silesia region obtained with the double functional median method (calculated with MBD).
Figure 6 presents PM10 concentration in the air forecast in calculated with Shang&Hyndman method (moving window equals 10) for five stations and for the Silesia region. Figure 7 presents five boxplots for the average sum of the differences between the observed curves and the curves forecasted with Shang&Hyndman method. Above a functional boxplot for the forecasts for Silesia region obtained with the Hyndman & Shang method
Our method seems to be more robust for functional outliers, what is not very surprising, as Shang and Hyndman make a forecasts basing on nonrobust generalized least squares method.
Predictor  Biel  Cze  Dab  Gli  Kat 

S&H  1699.585  1447.805  1699.585  881.2513  891.188 
Our forecasts  265.81  377.35  265.81  526.79  328.63 
5.2 Maximization of a social welfare
Generally speaking, it is known that air pollution has a negative impact on human health, however a form of this impact may take many different and very complex forms. Dangerous substances may interact one with another. A nature of impact may depend on age group and time of the day and night.
Our main aim is to maximize a "summarized" utility of a local community over a certain period, that symbolically may be written as (for details see Fleurbaey and Maniquet, 2011) :
(5.1) 
where is a number of a day, denotes social welfare related to PM10 emission reduction (positive and negative external effects valued in a fixed currency) and denotes a cost of the PM10 emission reduction valued in the fixed currency.
We assume that
and
It means that welfare related to PM10 is a function of an air quality (valued basing on evidenced costs of hospitalization due to lung diseases), medical expenses related to allergies (), an user friendless of a local environment (, a quality of a local information system providing information on air quality and health threats () and finally sociodemographic parameters of the community (). The cost related to the PM10 reduction relates to fixed costs including investments in new technologies (), variable costs involving effects tied with changes of weather causing lower or higher demand on heating energy (), political cost related to a transformation of popular heating systems basin for example on a coal into "clean" systems basing for example on nuclear energy, and costs related to quality of forecasting of the air pollution ().
In this paper we focus our attention on the last quantity measuring quality of forecasting of the air pollution in the selected region.
In our opinion it is reasonable to assume that the welfare associated with the pollution at considered town or region is linearly proportional to its population.
6 Conclusions
In the paper we have critically discussed an application of a model for hierarchical functional time series to studies of a day and night air pollution in Silesia region in Poland. We have focused our attention on
optimal estimation issues of the model. In this context we have compared our own original estimator with the best estimator known from the literature. We have assumed that an aggregate welfare of people living in the Silesia region is a function of among others a quality of the model estimation.
Our considerations clearly show usefulnesses of the HFTS methodology in a context of a local community welfare optimization. Shang and Hyndman approach provides elegant tools for HFTS modelling and forecasting in case of a relatively rich hierarchical structure and functional data without outliers. It is worth noticing, that our original double median HFTS predictor performs very well in a comparison to Shang and Hyndman predictor especially in terms of its computational complexity and robustness to functional outliers.
In our current research we concentrate on on the optimization issues related to specific forms of the formula (5.1) defining the welfare of a certain local community.
References

[1]
ArribasGil A., Romo J., (2014), Shape outlier detection and visualization for functional data: the outliergram, Biostatistics 15(4), 603–619.
 [2] Brown G.W., (1947), On SmallSample Estimation, Annals of Mathematical Statistics 18(4), 582585.
 [3] Ferraty F., and Vieu P., (2006), Nonparametric Functional Data Analysis: Theory and Practice, Springer.
 [4] Fliedner G., (2001), Hierarchical forecasting: issues and use guidelines, Industrial Management and Data Systems 101(1), 512.

[5]
Gijbels I., Nagy S., (2015), Consistency of nonintegrated depths for functional data, Journal of Multivariate Analysis 140, 259282.
 [6] Górecki T., Krzyśko M.K., Wołyński W., Wszak Ł, (2016), Selected statistical methods of data analysis for multivariate functional data, Statistical Papers, DOI10.1007/s0036201607578.
 [7] Horváth L., Kokoszka P., (2012), Inference for Functional Data with Applications, Springer, New York.
 [8] Hyndman R.J., Ahmed R.A., Athanasopoulos G., Shang H.L., (2011), Optimal combination forecasts for hierarchical time series, Computational Statistics & Data Analysis 55(9), 2579 2589.
 [9] Kahn K.B., (1998), Revisiting topdown versus bottomup forecasting, The Journal of Business Forecasting Methods & Systems 17(2), 1419.
 [10] Kohn R., (1982), When is an aggregate of a time series efficiently forecast by its past, Journal of Econometrics 18(3), 337349.
 [11] Kosiorowski D., (2014), Functional regression in shortterm prediction of economic time series, Statistics in Transition  new series 15(4), 611626.
 [12] Kosiorowski D., (2016), Dilemmas of robust analysis of economic data streams, Journal of Mathematical Sciences 1(2), 5972.
 [13] Kosiorowski D., Zawadzki Z. (2017) DepthProc: An R Package for Robust Exploration of Multidimensional Economic Phenomena, available at: arXiv:1408.4542v9.
 [14] Kosiorowski D., Rydlewski J.P., Snarska M., (2017a), Detecting a structural change in functional time series using local Wilcoxon statistic, Statistical Papers, DOI 10.1007/s003620170891y.
 [15] Kosiorowski D., Mielczarek D., Rydlewski J.P., Snarska M., (2017b), Generalized exponential smoothing in prediction of hierarchical time series, available at: arXiv:1612.02195v2.
 [16] Kosiorowski, D., Mielczarek, D., and Rydlewski, J. P., (2017c), Double functional median in robust prediction of hierarchical functional time series  an application to forecasting of the Internet service users behaviour, available at: arXiv:1710.02669v1.
 [17] Kosiorowski D., Rydlewski J.P., Zawadzki Z., (2017d), Functional outliers detection by the example of air quality monitoring, Statistical Review (in Polish, to appear).

[18]
Liu R., Parelius J.M., Singh K., (1999), Multivariate analysis by data depth: descriptive statistics, graphics and inference (with discussion), Annals of Statistics 27, 783858.
 [19] LópezPintado S., Romo J. ,(2007), Depthbased inference for functional data, Computational Statistics & Data Analysis 51(10), 4957–4968.
 [20] LópezPintado S., Romo J., (2009), On the concept of depth for functional data, Journal of the American Statistical Association 104, 718734.
 [21] LópezPintado S., Romo J., Torrente A., (2010), Robust depthbased tools for the analysis of gene expression data, Biostatistics 11(2), 254264.

[22]
Nagy S., Gijbels I., Omelka M., Hlubinka D., (2016), Integrated depth for functional data: Statistical properties and consistency, ESIAM Probability and Statistics 20, 95130.
 [23] NietoReyes A., Battey H., (2016), A Topologically Valid Definition of Depth for Functional Data, Statistical Science 31(1), 6179.
 [24] Pfanzagl J., (1970), On the Asymptotic Efficiency of Median Unbiased Estimates, The Annals of Mathematical Statistics, 41(5), 15001509.
 [25] Pfanzagl J., (1979), On optimal median unbiased estimators in the presence of nuisance parameters, The Annals of Statistics 7(1), 187193.
 [26] Ramsay J.O., Hooker G., Graves S., (2009), Functional Data Analysis with R and Matlab, Springer.
 [27] Shang H.L., Hyndman R.J., (2017), Grouped functional time series forecasting: an application to agespecific mortality rates, Journal of Computational and Graphical Statistics 26(2), 330343.
 [28] Shlifer E., Wolff R.W., (1979), Aggregation and proration in forecasting, Management Science 25(6), 594603.
 [29] Sun Y., Genton M.G., (2011), Functional Boxplots, Journal of Computational and Graphical Statistics 20(2), 316334.
 [30] Tabalerroni N., (2017), Robust Statistical Methods in Functional Data Analysis, Doctoral thesis and R package roahd, Politecnico di Milano.
 [31] Vinod H.D., Lopez de Lacalle J., (2009), Maximum entropy bootstrap for time series: the meboot R package, Journal of Statistical Software 29(5), 119.
 [32] Weale M., (1988), The reconciliation of values, volumes and prices in the national accounts, Journal of the Royal Statistical Society A 151(1), 211221.
 [33] Zuo Y., Serfling R., (2000), General notions of statistical depth function. Annals of Statistics 28(2), 461482.
 [34] Fleurbaey M., Maniquet F.,(2011), A Theory of Fairness and Social Welfare, Cambridge University Press.

[35]
http://powietrze.katowice.wios.gov.pl (access date: 24th March 2017).
WIOŚ underlines, that the data is gathered automatically and might be unverified