In Applied Statistics, the term *reconciliation* usually refers to adjusting several time series, corresponding to sectoral or geographically
disaggregated areas, to their aggregated total. Time series reconciliation is often
combined with another concept known as *temporal benchmarking*, which provides an adjustment of the estimated high-frequency series to known low-frequency
temporal series.

Even though the problem was first posed over seventy years ago (Deming & Stephan, 1940), new tools, and in particular new information technology tools, have sparked revived interest in the matter. A non-exhaustive list of old key references includes Di Fonzo and Filosa (1987), Cholette (1988), Chen and Dagum (1997), or Di Fonzo (1994, 2003).

Mention should be made of two highly innovative contributions. First, Di Fonzo and Marini (2005, 2011, 2015) used the so-called Denton’s movement preservation principle (Denton, 1971) and Causey and Trager’s (1981; available as an Appendix in Bozik & Otto [1988]) growth rates preservation principle to solve several reconciliation problems, including temporal or contemporaneous aggregation constraints, one or two-way time series systems or marginal benchmarking problems.

Second, the book by Dagum and Cholette (2006) analyses and solves benchmarking, calendarisation or reconciliation problems by using regression-based models. Said authors also derive solutions based on autoregressive integrated moving average (ARIMA) and structural time series models. Their approach includes distribution or interpolation problems. As in the Di Fonzo and Marini (2005) papers, the authors deal with one and two-way problems.

Present-day interest in studying the main topics concerning benchmarking and reconciliation for time series is evidenced by two recent publications: Firstly, the new version of the Quarterly National Accounts (QNA) Manual (International Monetary Fund [IMF], 2018) which dedicates chapter 6 to the systematic review of benchmarking methods that are relevant for compiling Quarterly National Accounts, although these methods may indeed prove helpful in addressing other problems, and secondly a special issue of Statistica Neerlandica (Chen, Di Fonzo, & Mushkudiani, 2018) devoted to the state-of-the-art of these topics. Several papers address specific aspects and tools (for example, Chen, Di Fonzo, Howells, & Marini 2018; Guerrero & Corona, 2018; or Bisio & Moauro, 2018, among others), whilst others offer a compilation of methods or procedures. For example, Quilis’ (2018) paper analyses different benchmarking procedures "in terms of practical feasibility, ease of use, and availability of dedicated software" (p. 448).

In this paper, the authors propose a method applicable to problems simultaneously involving both reconciliation and temporal benchmarking. The technique is herein applied to a Laspeyres-type volume index by employing a method derived from a proposal by Rojo and Sanz (2017), modified for use when cross-sectional restrictions employ weighted sums with time-varying weights. The solution is one-step, thereby simultaneously optimising benchmarking and reconciliation. The relevance of the problem is evident, for example, in the QNA Manual (IMF, 2018, pp. 194-195) which explores the possible inconsistences among aggregate QNA and their disaggregated estimates, with these inconsistencies deriving from the non-additivity of the Annual-Overlapping method. Said manual suggests reducing this inconvenience by presenting only percent measures of components’ contribution to the aggregated variable. The transversal non-additivity for the popular Annual-Overlapping method is a key challenge in cross-sectoral reconciliation for quarterly estimates of volume index.

The proposed method allows several indicators to be used, and does not require them to be approximations for the value to be estimated. It should also be pointed out that the stochastic nature of the model proposed by the authors enables the dispersion of the solution obtained to be estimated, thereby providing Bayesian confidence intervals, see Zellner (1971, pp. 27-28), for said solution, and also obtains the expression of the linear model which relates the indicators to the high-frequency series to be estimated.

Di Fonzo and Marini (2015) or Dagum and Cholette (2006) propose alternative methods, albeit focusing on the case in which (only) one indicator approximates the high-frequency series to be estimated. Cuevas et al. (2011, 2015) designed a method focused on benchmarking and reconciliation for National and Regional Accounts. These authors resolve the two aspects separately and, therefore, do not ensure that the final solution respects the set of restrictions. Other methodological differences concern the use of indicators. Although the Cuevas et al. (2015) method does allow for the use of multiple indicators, it initially combines them by means of dynamic factor analysis, such that in methodological terms it is a single indicator procedure like the previous ones.

In the following section, a detailed description is given for the assumptions and development of the proposed methodology, obtaining the explicit expression of the solution. The third section applies the proposed solution to obtain Quarterly Regional Accounts (QRA) for Spanish regions, with both the Annual National Accounts (ANA) and QNA for the whole of Spain being known. The method may obviously be applied to any country’s QNA and ANA. This example uses the same data as in the Cuevas et al. (2011) proposal, albeit over a longer time period. The final section summarises the most relevant findings. The work includes Supplementary Materials that contains colour tables and illustrations that are of collateral importance.

## Bayesian Benchmarking and Reconciliation in the Context of Time-Varying Aggregation

Let ${a}_{T}^{i},\text{\hspace{0.17em}}T=\mathrm{1,...,}N,\text{\hspace{0.17em}}i=\mathrm{1,...,}R$ be a non-stochastic annual variable, observed for $N$ consecutive years and for $R$ disaggregated entities (typically, the disaggregation is linked to a sectoral or geographic classification).

We also assume the annual variables resulting from the disaggregation over the classification,
${a}_{T}^{i}\text{,\hspace{0.17em}}T=\mathrm{1,...,}N$, to be known. Even though in the simplest applications the aggregation is achieved
through either the sum or the mean, we consider a more general scheme such that the
aggregate series is a *‘linear combination’* of the disaggregated ones, with the coefficients of the combination being time-varying,
however. Specifically, we suppose

##### 1

${a}_{T}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\displaystyle \sum _{i=1}^{R}{c}_{T}^{i}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{a}_{T}^{i}}\text{,\hspace{0.17em}}T=\mathrm{1,...,}N\text{\hspace{0.17em}}$The high or sub-annual frequency (usually monthly or quarterly)^{1} aggregated series is also assumed to be known,

where
$t,T$ refers to period
$t$ of year
$T$, and where the number of periods for a year will be denoted by
$m$. There may be an additional ‘*incomplete*’ year, in other words, with
$r$ sub-annual periods to be estimated,
${q}_{t,N+1}\text{,\hspace{0.17em}}t=\mathrm{1,...,}r$, for
$r<m$.

The aim is to estimate the sub-annual series for each disaggregated area,

including, if any, the high-frequency values relative to the incomplete year, ${q}_{t,N+1}^{i}\text{,\hspace{0.17em}}t=\mathrm{1,...,}r$, with $r<m$.

We assume that the cross-aggregation links at the sub-annual level are the same as are used for the annual series; that is to say

##### 2

${q}_{t,T}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\displaystyle \sum _{i=1}^{R}{c}_{T}^{i}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}{q}_{t,T}^{i}}\text{,\hspace{0.17em}}t=\mathrm{1,...,}m\text{,\hspace{0.17em}}T=\mathrm{1,...,}N\text{\hspace{0.17em}}$considering, if appropriate, the same relation for the incomplete year.

Finally, we should establish the temporal aggregation scheme for both levels of cross-aggregation. The arithmetic mean has been taken, although other classical schemes (sum, first or last values) are developed in an analogous fashion. Specifically, we assume

and

##### 3

$\text{\hspace{0.17em}}{a}_{T}^{i}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{1}{m}\text{\hspace{0.17em}}{\displaystyle \sum _{t=1}^{m}{q}_{t,T}^{i}}\text{,\hspace{0.17em}}i=\mathrm{1,...,}R,\text{\hspace{0.17em}}T=\mathrm{1,...,}N\text{\hspace{0.17em}}$All of the above classical schemes are consistent with the Equation 1 and Equation 2 cross-aggregation schemes, including the equality of coefficients ${c}_{T}^{i},\text{\hspace{0.17em}}T=\mathrm{1,...,}N\text{\hspace{0.17em}},i=\mathrm{1,...,}R$ for both expressions.

Before building the relevant distributions, we should stack the relevant series in matrix form. Specifically, we denote by ${q}^{i}*$ the column vector including the sub-annual chained series for the $i$-th individual area, i.e., for $i=\mathrm{1,...,}R$, ${q}^{i}*=({q}_{\mathrm{1,1}}^{i}\mathrm{,...,}{q}_{m,N}^{i})\text{'}$, if the last year is complete, and ${q}^{i}*=({q}_{\mathrm{1,1}}^{i}\mathrm{,...,}{q}_{r,N+1}^{i})\text{'}$ for incomplete last year schemes. We also denote by $q*$ the contemporaneously aggregated chain series $q*=({q}_{\mathrm{1,1}}^{}\mathrm{,...,}{q}_{m,N}^{})\text{'}$ or $q*=({q}_{\mathrm{1,1}}^{}\mathrm{,...,}{q}_{r,N+1}^{})\text{'}$ for both schemes.

In addition, we denote by $a*=({a}_{1}^{}\mathrm{,...,}{a}_{N}^{})\text{'}$ the column vector of annual chained series for the aggregated area, and by ${a}^{i}*$ the column vector of annual chained series for the $i$-th individual area, i.e., ${a}^{i}*=({a}_{1}^{i}\mathrm{,...,}{a}_{N}^{i})\text{'}\text{,\hspace{0.17em}}i=\mathrm{1,...,}R$. The cross-disaggregated series may then be stacked into the column vectors

Before continuing, one clarification concerning the dimensions of the matrices involved should be made; we denote by ${d}_{x}=NmR$ the column dimension of $x$, ${d}_{y}=NR$ the column dimension of $y$, and by ${d}_{q}=Nm$ the number of high-frequency periods. These dimension values are, respectively, ${d}_{x}=(Nm+r)R$, ${d}_{y}=NR$ and ${d}_{q}=Nm+r$ for incomplete year schemes. Obviously, ${d}_{x}=R{d}_{q}$. By using this notation, we can compactly write sub-annual data as

In sum, we have the annual and sub-annual series for the aggregated area, respectively, $a*$ and $q*$. We also have the annual series ${a}^{i}*$ for the different disaggregated areas, $i=\mathrm{1,...,}R$. The main aim is to estimate the sub-annual series for those disaggregated units, ${q}^{i}*\text{,\hspace{0.17em}}i=\mathrm{1,...,}R$, combined with the appropriate measures for the accuracy of the estimation.

As usual, Bayesian strategy initially states the prior distributions for the parameters and variables involved. These distributions include a certain number of non-random hyperparameters, whose values shall be stated by using allocation procedures. In addition, we then establish a behavioural linear model explaining the sub-annual variables as a function of one or more relevant indicators or approximated series. This linear model allows the likelihood function to be established and, consequently, the posterior distribution for the quarterly series and other relevant parameters to be derived.

### Joint Prior Distribution

Following Rojo and Sanz (2005, 2017), the authors propose a hierarchical Bayes normal-gamma approach. Specifically, we assume a normal conditional prior density for $x$, given the precision, $\tau $,

obeying restrictions (Equation 2) and (Equation 3), where
$\mu $ is the
$({d}_{x}\times 1)$ average vector for
$x$,
$\tau \text{\hspace{0.17em}}P$ is the
$({d}_{x}\times {d}_{x})$ precision matrix (the inverse of the variance-covariance matrix), and
$D={I}_{R}\otimes {D}^{\ast}$, with
${I}_{R}$ being the identity matrix of order
$R$ and
${D}^{\ast}$ the
$({d}_{q}-1)\times {d}_{q}$ first-difference matrix^{2} (with rank
${d}_{q}-1\text{\hspace{0.17em}}$). A greater degree of smoothness can be achieved by substituting
$D$ for
${D}_{2}={I}_{R}\otimes {D}_{2}^{\ast}$, with
${D}_{2}^{\ast}$ being the
$({d}_{q}-2)\times {d}_{q}$ matrix of rank
${d}_{q}-2$ that provides second order differences. Although only first order differences are
used in the formal derivation of the estimates, the results for second order differences
are obtained by simply replacing
$D\text{'}D$ with
${D}_{2}^{}\text{'}{D}_{2}$ in the following formulae.

We assume a gamma^{3} prior distribution for
$\tau $,

##### 4

$\pi \left(\tau \left|y\right.\right)\text{\hspace{0.17em}}\propto \text{\hspace{0.17em}}{\tau}^{\alpha -1}\text{\hspace{0.17em}}{e}^{-\tau \beta}\text{,\hspace{0.17em}}\tau >0\text{\hspace{0.17em}}$In sum, we obtain a normal-gamma prior joint distribution for $x$ and $\tau $,

##### 5

$\pi (x,\tau \left|y)\right.\propto \text{\hspace{0.17em}}{\tau}^{{d}_{x}/2+\alpha -1}\text{\hspace{0.17em}}\mathrm{exp}\left\{-\frac{\tau}{2}\left[2\beta +(x-\mu )\text{'}P(x-\mu )+x\text{'}D\text{'}Dx\right]\right\}\text{\hspace{0.17em}}$where $x$ obeys restrictions (Equation 2) and (Equation 3).

The prior distributions include a large number of parameters and, generally, a limited amount of data. These parameters are $\alpha $ and $\beta $ (from prior distribution for $\tau \text{\hspace{0.17em}}$), matrices $\overline{\delta}$ and ${P}_{\delta}\text{\hspace{0.17em}}$(involved in the prior for $\delta $), matrices $\mu $ and $P$ from the distribution for $x$ and, finally, matrix ${P}_{\epsilon}$ from the likelihood (see below).

Some authors suggest the use of Jeffreys’ rule to assign ‘*vague*’ or ‘*non-informative*’ priors. Examples include Box and Tiao (1973), Zellner (1971) or Gamerman (1997), among others. Other authors, Young (1996), Spiegelhalter et al. (1999), Broemeling (1985) or Rojo and Sanz (1999), proposed using ‘*almost*’ or ‘*approximately’* non-informative distributions as priors.

Non-informative priors are accurate for problems in which we could take large-sized
samples, allowing for a dominance of likelihood over priors. By contrast, for many
Macroeconomic and Psychological studies, or for numerous analyses related with Natural
Sciences, we have only a small amount of statistical information and, consequently,
prior distributions dominate likelihood. The authors have chosen the option of using
‘*informative’* priors, thus allocating reasonable values for the hyperparameters by using the information
included in the marginal prior for the relevant parameters and variables. Readers
can find the method in Rojo and Sanz (2005).

### Likelihood Function and Posterior Joint Distribution

In order to establish the likelihood function, we assume a linear model relating the disaggregated sub-annual series with ${k}_{i}$ indicators (maybe, thus, in different amounts for each disaggregated entity)

The model has a similar expression if the last year is incomplete.

We now write a matricial version for the proposed likelihood, which will provide more compact reasoning. We denote $\delta =({\delta}_{1}^{\text{'}}\mathrm{,...,}{\delta}_{R}^{\text{'}})\text{'}$ with ${\delta}_{i}=({\delta}_{\mathrm{1,}i}\mathrm{,...,}{\delta}_{{k}_{i},i})\text{'},\text{\hspace{0.17em}}i=\mathrm{1,...,}R$ the parameters of the likelihood model for each disaggregated area, and we include the values of the indicators in the matrix $Z=diag{\left\{{Z}_{i}\right\}}_{i=1}^{R}$, block diagonal matrix of size $\left({d}_{x}\times {\displaystyle \sum _{i=1}^{R}{k}_{i}}\right)$ whose elements are, for $i=\mathrm{1,...,}R$, ${Z}_{i}=\left({z}_{i\mathrm{,1}}\mathrm{,...,}{z}_{i,{k}_{i}}\right)$ being ${z}_{ij}=\left({z}_{1}^{ji}\mathrm{,...,}{z}_{{d}_{q}}^{ji}\right)\text{'}$.

Denoting the perturbations for the models as $\epsilon =({\epsilon}_{1}^{\text{'}}\mathrm{,...,}{\epsilon}_{R}^{\text{'}})\text{'}$, with ${\epsilon}_{i}=({\epsilon}_{{}_{1}}^{i}\mathrm{,...,}{\epsilon}_{{}_{{d}_{q}}}^{i})\text{'}$ the linear model could be written as

The likelihood is completely established by defining the prior for $\epsilon $ given $\tau $, $(\epsilon |\tau )$ as $N\left(0\text{,\hspace{0.17em}}{\tau}^{-1}\text{\hspace{0.17em}}{P}_{\epsilon}^{-1}\right)$, assuming that $Cov({\epsilon}_{i},{\epsilon}_{j})=\mathrm{0,}\text{\hspace{0.17em}}i\ne j$ and $VCov({\epsilon}_{i})={P}_{\epsilon}^{i}$, and being the hyperparameter ${P}_{\epsilon}$ a block diagonal matrix, ${P}_{\epsilon}\text{\hspace{0.17em}}=diag{\left\{{P}_{\epsilon}^{i}\right\}}_{i=1}^{R}$.

Hence, the likelihood function is given by

##### 6

$L(x,\delta ,\tau |y,Z)\text{\hspace{0.17em}}\propto \text{\hspace{0.17em}}{\tau}^{{d}_{x}/2}\text{\hspace{0.17em}}\mathrm{exp}\left\{-\frac{\tau}{2}(x-Z\delta )\text{'}{P}_{\epsilon}(x-Z\delta )\right\}$being $y$ the annual data vector for the disaggregated areas previously defined.

We assume a normal prior distribution for $\delta $, given $\tau $,

##### 7

$\pi (\delta |\tau ,y)\text{\hspace{0.17em}}\propto \text{\hspace{0.17em}}{\tau}^{k/2}\text{\hspace{0.17em}}\mathrm{exp}\left\{-\frac{\tau}{2}(\delta -\overline{\delta})\text{'}{P}_{\delta}(\delta -\overline{\delta})\right\}$with $k={\displaystyle \sum _{i=1}^{R}{k}_{i}}$, and ${P}_{\delta}$ being the precision matrix, i.e., the conditional prior for $\delta $ is a ${N}_{k}\left(\overline{\delta}\text{,\hspace{0.17em}}{\tau}^{-1}\text{\hspace{0.17em}}{P}_{\delta}^{-1}\right)$.

Using Equation 5, Equation 6 and Equation 7, and assuming prior independence between $\delta $ and the remaining parameters, except $\tau $, the posterior distribution can be expressed as

##### 8

$$\begin{array}{l}p(x,\delta ,\tau |y,Z)\text{\hspace{0.17em}}\propto \text{\hspace{0.17em}}{\tau}^{{d}_{x}+\alpha -1+k/2}\text{\hspace{0.17em}}\cdot \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdot \mathrm{exp}\left\{-\frac{\tau}{2}\right.\left[2\beta +(x-\mu )\text{'}P(x-\mu )+x\text{'}D\text{'}Dx+\right.\left.\left.(\delta -\overline{\delta})\text{'}{P}_{\delta}(\delta -\overline{\delta})+(x-Z\delta )\text{'}{P}_{\epsilon}(x-Z\delta )\right]\right\}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\end{array}$$under restrictions (Equation 2) and (Equation 3).

A well-known result (Zellner, 1971, p. 30) states that for a quadratic loss function

with $C$ being a positive definite matrix, the minimum of the quadratic risk function (the average of the quadratic loss function with the posterior distribution) is achieved by taking for $\left(\tau ,x,\delta \right)$ the mean of that posterior joint distribution. The solution, therefore, involves obtaining said posterior means.

Note, however, that the restricted distribution for
$x$ is a degenerate one. Obtaining the posterior average for
$x$ may be achieved by including the temporal and contemporaneous restrictions in the
posterior distribution. This substitution removes several components of
$x$. We then obtain the posterior mean for the ‘*active*’ parameters and variables and, taking into account the linearity of the restrictions,
we derive the outcome for the remaining parameters.

Specifically, our aim is to write $x=W{x}^{r}+{y}_{v}$, where $W$ and ${y}_{v}$ are non-random matrix and vector, respectively, and ${x}^{r}$ will group the components that are not subject to restrictions for the sub-annual series. The temporal constraints (Equation 1) lead to one sub-annual period being excluded for each year (specifically, we have excluded the last sub-annual period for each year). Furthermore, the transversal aggregation constraint (Equation 2) implies the linear dependence of one disaggregated sub-annual series. Thus, the $R$-th series has also been excluded.

Denoting by
${q}^{i,r}$(
$r$ for ‘*restricted’*) the column vector whose components are the ‘*independent*’ values for the chained series corresponding to the
$i$-th disaggregated area (
${q}^{i,r}=({q}_{\mathrm{1,1}}^{i}\mathrm{,...,}{q}_{m-\mathrm{1,1}}^{i}\text{\hspace{0.17em}}\mathrm{,...,}\text{\hspace{0.17em}}{q}_{\mathrm{1,}N}^{i}\mathrm{,...,}{q}_{m-\mathrm{1,}N}^{i})\text{'}$ or
${q}^{i,r}=({q}_{\mathrm{1,1}}^{i}\mathrm{,...,}{q}_{m-\mathrm{1,1}}^{i}\text{\hspace{0.17em}}\mathrm{,...,}\text{\hspace{0.17em}}{q}_{\mathrm{1,}N}^{i}\mathrm{,...,}{q}_{m-\mathrm{1,}N}^{i},{q}_{\mathrm{1,}N+1}^{i}\mathrm{,...,}{q}_{r,N+1}^{i})\text{'}$, depending on the presence of an incomplete last year), these independent values
can be grouped into the column vector
${x}^{r}=({q}^{\mathrm{1,}r}\text{\hspace{0.17em}}\text{'}\mathrm{,...,}{q}^{R-\mathrm{1,}r}\text{\hspace{0.17em}}\text{'})\text{'}$. Note that the dimension of
${x}^{r}$ is equal to
${d}_{{x}_{r}}={d}_{x}-{d}_{q}-N(R-1)$.

Note S1 in Supplementary Materials provides the linear relation linking $x$ and ${x}^{r}$

##### 9

$x\text{\hspace{0.17em}}=\text{\hspace{0.17em}}W\xb7{x}^{r}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{y}_{v}$Now, integrating density (Equation 8), we obtain the posterior marginal distributions for $x$ and for $\delta $ and, replacing restriction Equation 9 provides us with the posterior restricted distributions for both ${x}^{r}$ and $\delta $,

##### 10

$p({x}_{r}|y,Z)\text{\hspace{0.17em}}\propto \text{\hspace{0.17em}}\text{\hspace{0.17em}}{\left[1+\left({x}_{r}-{a}_{r}\right)\text{'}\frac{W\text{'}MW}{{b}_{r}}\left({x}_{r}-{a}_{r}\right)\right]}^{-({d}_{x}+\alpha )}$being $M=P+D\text{'}D+{P}_{\epsilon}-{P}_{\epsilon}Z{P}_{\delta \epsilon}^{-1}Z\text{'}{P}_{\epsilon}$ a square matrix, with ${P}_{\delta \epsilon}={P}_{\delta}+Z\text{'}{P}_{\epsilon}Z$, ${b}_{r}=b+\left({M}^{-1}u-{y}_{v}\right)\text{'}M\left[{M}^{-1}-W{\left(W\text{'}MW\right)}^{-1}W\text{'}\right]M\left({M}^{-1}u-{y}_{v}\right)$, for $b=s-u\text{'}{M}^{-1}u$ and $s=2\beta +\mu \text{'}P\mu +\overline{\delta}\text{'}{P}_{\delta}\overline{\delta}-\overline{\delta}\text{'}{P}_{\delta}{P}_{\delta \epsilon}^{-1}{P}_{\delta}\overline{\delta}$ scalars, and $u={P}_{\epsilon}Z{P}_{\delta \epsilon}^{-1}{P}_{\delta}\overline{\delta}+P\mu $ a vector, being ${a}_{r}={\left(W\text{'}MW\right)}^{-1}W\text{'}M\left({M}^{-1}u-{y}_{v}\right)$.

Posterior density (Equation 10) is a Student multivariate (henceforth, MS-t)^{4} with
${\upsilon}_{{x}_{r}}\text{\hspace{0.17em}}={d}_{q}(1+R)+2\alpha +N(R-1)$ degrees of freedom. The scale matrix is
$\text{\hspace{0.17em}}\frac{{\upsilon}_{{x}_{r}}}{{b}_{r}}W\text{'}MW$ and the position vector
${a}_{r}$. Thus, the posterior variance-covariance matrix for
${x}_{r}$ could be written as
$V({x}_{r})=\frac{{\upsilon}_{{x}_{r}}}{{\upsilon}_{{x}_{r}}-2}\text{\hspace{0.17em}}\frac{{b}_{r}}{{\upsilon}_{{x}_{r}}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{(W\text{'}MW\text{\hspace{0.17em}})}^{-1}$.

The posterior first and second order moments for $x$ are then

##### 11

$E(x|y,Z)=W{a}_{r}+{y}_{v}$and

##### 12

$V(x)=\frac{{b}_{r}}{{\upsilon}_{{x}_{r}}-2}W{(W\text{'}MW)}^{-1}W\text{'}$Posterior density for $\delta $ is obtained in a similar manner. We obtain

##### 13

$$\pi (\delta )\propto {\left[{b}_{\delta}+(\delta -{a}_{\delta})\text{'}{M}_{\delta}(\delta -{a}_{\delta})\right]}^{\text{\hspace{0.17em}}-\text{\hspace{0.17em}}\frac{k+(2{d}_{x}-{d}_{{x}_{r}}+2\alpha )}{2}}$$with ${M}_{\delta}={P}_{\delta}+Z\text{'}{P}_{\epsilon}({P}_{\epsilon}^{-1}-{Q}_{r}){P}_{\epsilon}Z$ a square matrix with dimension $k$ for ${Q}_{r}=W{(W\text{'}{P}_{r}W)}^{-1}W\text{'}$, a singular square matrix with dimension ${d}_{x}$ and rank ${d}_{{x}_{r}}$and being ${P}_{r}=P+D\text{'}D+{P}_{\epsilon}$, and with

and being ${a}_{\delta}$ the column vector ${a}_{\delta}={M}_{\delta}^{-1}\left[{P}_{\delta}\overline{\delta}+Z\text{'}{P}_{\epsilon}\left({y}_{v}+{Q}_{r}(P\mu -{P}_{r}{y}_{v})\right)\right]$.

The posterior distribution (Equation 13) for $\delta $ is, thus, an MS-t with $2{d}_{x}-{d}_{{x}_{r}}+2\alpha $ degrees of freedom, position ${a}_{\delta}$ and scale matrix $(2{d}_{x}-{d}_{{x}_{r}}+2\alpha )\frac{{M}_{\delta}}{{b}_{\delta}}$. The variance-covariance matrix is then equal to $V(\delta )=\frac{{b}_{\delta}}{2{d}_{x}-2{d}_{{x}_{r}}+2\alpha -2}{M}_{\delta}^{-1}$.

## An Example From the Quarterly Regional Accounts

We now present an example, which illustrates the above method. The resulting estimated
variables are the quarterly regional chained volume series for Spanish regional Gross
Domestic Product (GDP)^{5}.

Temporal restrictions impose consistency among the regional estimations (quarterly regional chained volume series) and annual chained volume series provided by the annual regional accounts (ARA) from official Spanish regional statistics.

We also force the transversal consistency of these quarterly regional series with the Spanish national chained series provided by the QNA, a consistency which states that the QNA chained series is a weighted sum of the regional chained volume indices.

Many national statistics institutes have used annual chain-linking series for ANA and corresponding quarterly series for QNA. Specifically, the US Bureau of Economic Analysis (BEA) has used quarterly chain-linking volume series since 1996. In the European Union (EU), a European task force was set up in 2007, co-chaired by Eurostat and the European Central Bank. The growing popularity of chain-linking series for both QNA and ANA has led to the need for efficient tools in the reconciliation and benchmarking of quarterly chain-linking series. Cuevas et al. (2011, 2015), for example, proposed a two-step-method for use with the annual overlap derivation of chain-linked QNA, a very popular technique frequently used in the EU (see Eiglsperger, 2008 for a detailed listing of its dissemination).

Chen, Di Fonzo, Howells, and Marini (2018) developed approaches for reconciling annual (preliminary) estimates of US national accounts aggregates subject to quinquennial benchmarks available from detailed input-output tables. Furthermore, in an update of the 2001 version of the IMF QNA manual (Shrestha [2013], a document prepared by M. Marini and Th. Alexander), the authors suggest that "the compilation of consistent quarterly estimates satisfying both low frequency benchmarks and accounting identities at the quarterly level has become more and more challenging for compilers" (Shrestha, 2013, p. 7). Although the UN’s System of National Accounts (SNA) recommends use of the Fisher-type index, the authors of the manual recognise that the Laspeyres volume indices are an acceptable alternative in national accounts.

More specifically, the Spanish National Statistics Institute (INE) provides the annual
GDP (chain-linked volumes, reference year 2010^{6}) for the 17 autonomous regions and for the two autonomous cities (hereinafter, 19
regions) for the nineteen-year period 2000-2018. The INE also estimates total quarterly
GDP (from QNA), both at market prices (in euros) and by the volume-chained series
(the annual overlap method is used here). In both cases, the raw series as well as
the seasonally and working-day adjusted series (SA) are presented. At the time of
writing this paper, the quarterly series was composed of 76 quarters^{7}, from the first quarter of 2000 to the fourth quarter of 2018.

As mentioned before, the authors’ aim is to estimate the 19-quarterly regional chained series, all being consistent with the annual regional chained series and with the total quarterly one. We only estimate the SA regional series, with the estimation for raw series following a similar development. We apply the procedure obtained in the second section for the period from quarter 2000:1 to 2018:4.

Cuevas et al. (2011) solve this problem by using a multivariate two-step extension of the Denton (1971) method. They use a procedure proposed by Di Fonzo and Marini (2005). The former authors use total employment (regional social security contributors [SSC]) as high-frequency indicators, among others, due to its close linkage to real output. Although the widespread debate concerning the extent of cyclical synchrony between employment and output is well known, we do not lead or lag the SSC series in order to replicate their data selection. As already pointed out, Cuevas et al. (2011) use a two-stage procedure that is not simultaneous, unlike the one proposed in this work.

Taking into account that the estimated quarterly regional series will be seasonally adjusted, the SSC series have first been seasonally adjusted using the X12 method, as implemented in Eviews (Version 6) software. The procedure proposed by the authors was implemented in MATLAB (Version R2012b).

Before the method can be applied, an extra adjustment is needed. This is due to the
lack of consistency among the whole set of annual regional chained-series and the
total Spanish one, resulting from the existence of ‘*extra-regio*’ territories (extra-regional economic activities such as embassies, military or scientific
bases or resource deposits in international waters, among others).

To sum up, we have regional series concerning the annual GDP for the 19 regions for the period 2000-2018, both at current prices and in terms of volume. We have also estimated Spanish quarterly GDP (without extra-regio), both raw series and SA series, also at current prices and in chained-linking terms. Finally, we know the raw SSC series for each region, and have previously derived the corresponding SA series. As pointed out earlier, our objective is to estimate the Quarterly Regional chained volume series for Spanish regional GDP.

We then apply the method obtained in the second section, whose notations we now describe.

The INE provides regional annual GDP at market prices, in euros, ${v}_{T}^{i}\text{,\hspace{0.17em}}T=\mathrm{1,...,}19\text{,\hspace{0.17em}}i=\mathrm{1,...,19}$, for the 17 autonomous regions and the two autonomous cities ( $R=19$ disaggregated areas) and for $N=19$ years (from 2000 to 2018). National annual GDP at market prices ${v}_{T}^{}\text{,\hspace{0.17em}}T=\mathrm{1,...,}19$ is obtained by aggregation. It also provides the annual volume chain series, ${a}_{T}\text{,\hspace{0.17em}}{a}_{T}^{i}\text{,\hspace{0.17em}}T=\mathrm{1,...,}19\text{,\hspace{0.17em}}i=\mathrm{1,...,19}$ at the national and regional level, respectively. Furthermore, the QNA provide the national chained volume series, ${q}_{t,T}\text{,\hspace{0.17em}}T=\mathrm{1,...,19}\text{,\hspace{0.17em}}t=\mathrm{1,...,4}$, seasonally and working-day adjusted. Finally, the Labour, Migration & Social Security Ministry provides the indicator used, ${Z}_{i}={\left\{{z}_{t,T}^{i}\right\}}_{T=\mathrm{1,...,19,}t=\mathrm{1,...,4}}\text{,\hspace{0.17em}}i=\mathrm{1,...,19}$, the social security contributors, SSC, seasonally corrected by the authors.

The ‘*annual overlap’* method employed by the INE establishes the links between high frequency (quarterly)
and low-frequency (annual) volume indices by using arithmetic means, both for national
and for regional volume index. We thus impose temporal benchmarking

Furthermore, the cross-restrictions are

##### 14

$${a}_{T}={a}_{T-1}\cdot {\displaystyle \sum _{i=1}^{19}{\omega}_{T-1}^{i}}\cdot \frac{{a}_{T}^{i}}{{a}_{T-1}^{i}}\text{,\hspace{0.17em}}T=\mathrm{1,...,19}$$at the annual level and

##### 15

$${q}_{t,T}={a}_{T-1}\cdot {\displaystyle \sum _{i=1}^{19}{\omega}_{T-1}^{i}\left(\frac{{q}_{t,T}^{i}}{{a}_{T-1}^{i}}\right)}\text{,\hspace{0.17em}}t=\mathrm{1,...,4}\text{,\hspace{0.17em}}T=\mathrm{1,...,19}$$at the quarterly one, being

Readers may easily note that the above Equation 14 and Equation 15 can be written as

with ${c}_{T}^{i}={a}_{T-1}\text{\hspace{0.17em}}{\omega}_{T-1}^{i}\text{\hspace{0.17em}}/\text{\hspace{0.17em}}{a}_{T-1}^{i}$, $T=\mathrm{1,...,19,}\text{}i=\mathrm{1,...,19}$. We are, therefore, in the context foreseen in section 2, such that the method proposed allows us to estimate the chained volume series ${q}_{t,T}^{i}\text{,\hspace{0.17em}}T=\mathrm{1,...,19}\text{,\hspace{0.17em}}t=\mathrm{1,...,4}\text{,\hspace{0.17em}}i=\mathrm{1,...,19}$, grouped in vector $x$, and to obtain its variance-covariance matrix.

Some representative tables and charts for the regional chain-linking SA series are shown in the Supplementary Materials. The authors have also obtained the results for the regional raw series, although these are not displayed in the work for reasons of length. Readers may observe the similarities and differences between our tables and charts and those obtained by Cuevas et al. (2011) for another time interval. These are shown in said reference.

Table S1 in the Supplementary Materials presents the Pearson correlations between SSC and estimated series, and also between
the growth rates (annual grow rates,
${t}_{(\mathrm{1,1})}=({y}_{t}-{y}_{t-1})/{y}_{t-1}$ and quarterly growth rates,
${t}_{(\mathrm{4,1})}=({y}_{t}-{y}_{t-4})/{y}_{t-4}$) for both series. Broadly speaking, the highest values are obtained for annual growth
rates, except for Ceuta and Melilla^{8} (the two autonomous cities). One surprising result is the low correlation in levels
for Galicia, the Canary Islands, and the Comunidad Valenciana.

Tables S2 to S4 in the Supplementary Materials provide an overview of the results obtained for the Spanish regions, showing in particular regional growth behaviour during the recent Great Recession.

For its part, Figures S1 and S2 offer a graphic overview of said regional growth path.

Table S5 in the Supplementary Materials presents the comparison between the cyclical signal of each region and the turning points for the aggregate Spanish reference by using ratios, as defined in Abad and Quilis (2004). Only the turning points are taken into account, and no specific attention is paid to the numerical values for the two series.

The first and second column show the so-called ‘conformity ratio’, comparing the paired^{9} turning points of each regional cyclical signal with the individual turning points.
R_{x} thus compares such paired turning points as a percentage of the turning points for
regional cycles, and R_{y} compares them as a percentage of national ones. The conformity ratio varies between
0% and 100%, showing the extent to which the paired turning points reflect the overall
cyclical signal of the region.

Readers may note that, with some exceptions, the agreement between regional and national cycles shown by indices Rx and Ry is relevant, with the exceptions corresponding to the Spanish regions of Ceuta and Castilla La Mancha.

The third column shows the global median delay (GMD) between regional and national cycles. Thus, the series are classified as coincident, lagged or leading with respect to the national cycle, respectively, for small (we take between -1 and +1), positive or negative GMD values. It should be noted that the Balearic Islands lead the national economy by two and a half quarters, and that Cantabria lags the national economy by two quarters.

The last column shows an index of cyclical coincidence^{10} between regional and national series. Values near 1 suggest that regional and national
series are highly procyclical, being highly countercyclical if the index approaches
the value -1. For values near zero, we conclude that the regional series is non-classifiable
compared to the national one. It should be noted that all of the values are positive,
but that only Andalusia, Cataluña, the Comunidad Valenciana, Madrid and the Basque
Country display high values.

We also performed the estimation with incomplete years. Now $T=18$ (between 2000 and 2017). The regions involved do not change, and we take $r=3$; in other words, we estimate the chained volume series for the three first quarters of 2018. The estimates are shown only for one region (Andalucía) in Figure S3 in the Supplementary Materials. Table S6 shows the Bayesian confidence intervals at 95% for this estimation (see Supplementary Materials).

## Conclusions

In this article, the authors propose a method to obtain explicit solutions for simultaneous benchmarking and reconciliation problems for a system of time series when the cross-restrictions use time-varying coefficients. The method provides explicit solutions to the estimation problem and deals with concurrently solving temporal restrictions (benchmarking annual and sub-annual frequency series) and contemporaneous ones (reconciliation among disaggregated and aggregated sub-annual frequency series).

The Bayesian model involved belongs to the frequently used normal-gamma family and minimizes a risk function derived from a quadratic loss function. In addition, the design of the method allows users to include one or several performance indicators through the likelihood model, and to estimate quarterly values for incomplete years.

The stochastic nature of the proposed model allows Bayesian confidence intervals to be obtained for each of the values of the high-frequency series estimated. These intervals are particularly interesting when estimating incomplete years.

Comparisons with alternative methods are not, broadly speaking, feasible since the methods have a different statistical base (mathematical methods compared to statistical methods) and are not nested models (none of them is a generalisation of the other). Nevertheless, certain differences may be pointed out.

Compared to the proposal put forward by Cuevas et al. (2011, p. 6, 2015, p. 631) it is worth highlighting that our procedure allows for the use of several indicators when estimating the sub-annual disaggregate series. Said indicators need not even be the same for all of them. In contrast, the proposal of Cuevas et al., previously synthesises the indicators in a single indicator through dynamic factor analysis. In addition, these authors’ procedure initially resolves temporal benchmarking and, subsequently, reconciliation of the disaggregated areas with the total aggregate. It is therefore (in any case) a sequentially optimal procedure, whereas the one proposed in this paper is globally optimal.

As regards the procedure of Di Fonzo and Marini (2011, 2015), we highlight that the design of their model implies using a single indicator. In
fact, said authors follow the original idea of Denton (1971), such that the problem is one of adjustment through a single indicator, which is
an approximation of the target series (for example, in Di Fonzo & Marini, 2011), they state that: "This procedure performs the constrained optimization of an objective
function according to which the proportionate difference between the benchmarked and
*the original series ... must be as constant as possible through time*" (emphasis ours; p. 148). In the words of IMF (2018), "the objective is to combine the quarterly movements of the indicator with the annual
levels of the ANA variables" (p. 87). Strictly speaking, therefore, it involves not
having an indicator but rather an approximation to the sub-annual series to be estimated.

Since it is an approximation, Di Fonzo and Marini (2011) suggest that one of the advantages of their method involves the similarity of the series estimated, whether either one or more, to the approximation (indicator) used. This apparent advantage proves to be a drawback when the indicator used as an approximation is too volatile, as tends to happen for small disaggregated areas. As seen in the second section, the proposed inclusion of an additional factor in the prior distribution seeks to correct, at least partially, that volatility, should it occur. The authors have simulated examples at different levels of volatility for the indicators, and have shown the ability of the Bayesian proposal to obtain smooth estimates. The comparison is carried out only for temporal benchmarking procedure and not for the reconciliation, given that Di Fonzo and Marini did not provide solutions for the reconciliation when time-varying links are present.

An example related to Spanish quarterly accounts, combining national and regional volume series, is presented, and evidences the method’s feasibility and appropriateness. In addition, individual regional behaviour during the recent twin economic crisis is analysed using the estimated quarterly volume series.

As previously pointed out (and illustrated in Figure S3 and Table S6 in the Supplementary Materials), the proposed method allows, as in those of Di Fonzo and Marini (2011) and Dagum and Cholette (2006), the target series to be estimated when the final year is incomplete. The additional advantage of the proposal put forward is that it enables Bayesian confidence bands to be built based on the posterior distribution of the estimated series.