Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
58,170 characters · 12 sections · 61 citation commands
A robust score-driven filter for multivariate time series
\def\spacingset#1{ {#1}} \spacingset{1}
\if00 \fi
\if10 {
} \fi
{\it Keywords:} Robust filtering, Multivariate models, Score-driven models, Homescan data.
\spacingset{1.2}
The analysis of multivariate time series has a long history, due to the empirical evidence, from most research fields, that time series resulting from complex phenomena do not only depend on their own past, but also on the history of other variables. For this reason, from Hannan1970, the literature on multivariate time series has grown very fast. The leading example is the dynamic representation of the conditional mean of a vector process which gives rise to vector autoregressive processes, see Hamilton1994 and Lutkepohl2007.
Following the taxonomy proposed in Cox1981, two main classes of models can be considered when analysing dynamic phenomena: parameter-driven and observation-driven models. The class of parameter driven model is a broad class, which involves unobserved component models and state space models (Harvey1989, West_Harrison1997). Within this framework, parameters are allowed to vary over time as dynamic processes driven by idiosyncratic innovations. Hence, likelihood functions are analytically tractable only in specific cases, notably linear Gaussian models, where inference can be handled by the Kalman filter. On the other hand, parameter-driven models are very sensitive to small deviations from the distributional assumptions. In addition, the Gaussian assumption often turns out to be restrictive, and flexible specifications may be more appropriate. Thus, a fast growing field of research is dealing with nonlinear or non-Gaussian state-space models, resting on computer intensive simulation methods like the particle filter discussed in durbkoopman2012. Although these methods provide extremely powerful instruments for estimating nonlinear and/or non-Gaussian models, they can be computationally demanding. Furthermore, it may be difficult to derive the statistical properties of the implied estimators, due to the complexity of the joint likelihood function.
In contrast, in observation-driven models, the dynamics of time varying parameters depend on deterministic functions of lagged variables. This enables a stochastic evolution of the parameters which become predictable given the past observations. Koopman_Lucas_Scharth2016 assess the performances and optimality properties of the two classes of models, in terms of their predictive likelihood. The main advantage of observation-driven models is that the likelihood function is available in closed form, even in nonlinear and/or non-Gaussian cases. Thus, the asymptotic analysis of the estimators becomes feasible and computational costs are reduced drastically.
Within the class of observation-driven models, score-driven models are a valid option for modeling time series that do not fall in the category of linear Gaussian processes. Examples have been proposed in the context of volatility estimation and originally referred to as generalised autoregressive score (GAS) models, CrealKoopLucas2013, and as dynamic conditional score (DCS) models, Harvey2013. The key feature of these models is that the dynamics of time-varying parameters are driven by the score of the conditional likelihood, which needs not necessarily to be Gaussian but can be heavy tailed. For example, it may follow a Student's t distribution as in HarveyLuati2014 and Linton2020, an exponential generalized beta distribution, as in Caivano_Harvey_Luati2016, a binomial distribution as in the vaccine example by hansen.2019, or represented by a mixture, see Lucas_Sch_Sch_2019. The optimality of the score as a driving force for time varying parameters in observation-driven models is discussed in Blasques_Koopman_Lucas2015. According to which conditional distribution is adopted, specific situations may be conveniently handled due to the properties of the score. As an example for the univariate case, if a heavy-tailed distribution is specified, namely Student's t, the resulting score-driven model yields a simple and natural model-based signal extraction filter which is robust to extreme observations, without any external interventions or diagnostics, like dummy variables or outlier detection, see HarveyLuati2014.
In score-driven models, as well as in all observation-driven models, the time varying parameters are updated by filtering procedures, i.e. weighted sums of functions of past observations, given some initial conditions that can be fixed or estimated along with the static parameters. A robust filtering procedure should assign less weight to extreme observations in order to prevent biased inference of the signal and the parameters. In particular, the work of Calvet_Czellar_Ronchetti2015 provides a remarkable application of robust methods when dealing with contaminated observations. The authors show that a substantial efficiency gain can be achieved by huberizing the derivative of the $\log$-observation density. As we show in the present study, the same holds if one considers an alternative robustification method, based on the specification of a conditional multivariate Student's t distribution. A similar approach can be found in Prucha_Kelejian1984 and Fiorentini_Sentana_Calzolari2003, where the multivariate Student's t distribution provides a valid alternative to relax the normality assumption. In the context of score-driven models, creal2014 mention the relevance of modeling high-frequency data with outliers and heavy tails by means of the multivariate Student's $t$ distribution.
In this paper, we develop a score-driven filter for the time-varying location of a multivariate Student's t distribution. The specification is similar to the multivariate model for the location addressed in Harvey2013 and has some traits in common with the quasi-vector autoregressive model by Blazsekwp, though our perspective is more focused on the aspects of the filter and its stochastic properties. A spatial extension of the model developed in this paper is considered by Gasperoni2021. We envisage three main contributions to the existing literature.
The first contribution is the derivation of the probabilistic theory behind the multivariate dynamic score-driven filter for conditional Student's t distributions, including the conditions of stationarity, ergodicity and invertibility, in a similar spirit of Comte2003 and Hafner2009 for the multivariate conditional variance models by baba1990 and Engle1995. These results provide the basis for further generalisations, such as, for instance, the spatial model by Gasperoni2021. As the conditional likelihood is available in close form, we estimate the static parameters with the method of maximum likelihood and prove strong consistency and asymptotic normality of the estimators. It is noteworthy to remark that when the degrees of freedom of the Student's $t$ distribution tend to infinity, we recover a linear Gaussian state-space model.
The second contribution is the development of an estimation scheme grounded on Fisher's scoring method, based on closed-form analytic expressions, which can be directly implemented into any statistical or matrix-friendly software.
The third contribution of the paper is an innovative application, dealing with estimation of regional consumer prices based on home scanner data. The use of scanner data to compute official consumer price indices (CPIs) is gaining popularity, because of their timeliness and a high level of product and geographical detail Feenstra2003. However, they also suffer from a variety of shortcomings, which make time series of scanner data prices (SDPs) potentially very noisy, especially when they are estimated for population sub-groups, or at the regional level Silver1995. There is extensive research and a lively debate on the issues related to the computation and use of scanner data based CPIs. In a dedicated session of the 2019 meeting of the the Ottawa Group on Price Indices, it has been suggested\footnote{See Jens Mehroff presentation at \url{https://eventos.fgv.br/sites/eventos.fgv.br/files/arquivos/u161/towards_a_new_paradigm_for_scanner_data_price_indices_0.pdf}} to adopt model-based filtering techniques to extract the signal from scanner-based time series of price data. These filtered estimates lose the classical price index formula interpretation, but are expected to deliver the same information content with a better signal-to-noise ratio. We show that our robust multivariate model, applied to SDPs, provides information on the dynamics of the time series and on their interrelations without being affected from outlying observations, which are naturally downweighted in the updating mechanism.
The paper is organised as follows. In Section (ref) the filter is specified. Section (ref) deals with the stochastic properties of the filter, while in Section (ref), likelihood inference is discussed. The empirical analysis is reported in section (ref). Some concluding remarks are drawn in Section (ref). The proofs of the results stated in the paper are collected in Appendix A. Online supplementary materials contain the details of a Monte Carlo study designed to assess the finite sample properties of the estimators, the relevant quantities for the implementation of the Fisher scoring algorithm as well as the proofs of some auxiliary Lemmata.
Let us consider a $\mathbb{R}^N$-vector of stochastic processes $\{ \boldsymbol{y}_t\}_{t\in\mathbb{Z}}$, $N \geq 1$, and let $\mathcal{F}_{t-1} = \sigma\{ \boldsymbol{y}_{t-1}, \boldsymbol{y}_{t-2}, $ $\boldsymbol{y}_{t-3}, \dots \}$ be its filtration at time $t-1$. The following stochastic representation of $\boldsymbol{y}_t$ is considered,
where $\boldsymbol{\mu}_t$ is a time varying location vector of $\mathbb{R}^N$, $\boldsymbol{\Omega}$ is a $N\times N$ scale matrix that we assume to be static and $\boldsymbol{\epsilon}_t \sim \boldsymbol{t}_\nu (\boldsymbol{0}_N, \boldsymbol{I}_N)$ is an independent identically distributed (IID) multivariate standard t-variate. With $\boldsymbol{0}_N$ we denote the null vector of $\mathbb{R}^N$ and with $\boldsymbol{I}_N$ the $N \times N$ identity matrix.
Our interest is in recovering $\boldsymbol{\mu}_t$ based on a set of observed time series from $\boldsymbol{y}_t$, for $t=1, \dots, T$, where $T \in \mathbb{N}$. With no distributional assumptions on the dynamics of $\boldsymbol{\mu}_t$, a filter can be specified,
that is a stochastic recurrence equation (SRE), where $\boldsymbol{\theta}\in \boldsymbol{\Theta} \subset \mathbb{R}^p$ is a vector of unknown static parameters, $\boldsymbol{\mu}_{t|t-1}$ is a $\mathbb{R}^N$-random vector that takes values in $\boldsymbol{\mathcal{M}} \subset \mathbb{R}^N$ and $\phi: \boldsymbol{\mathcal{M}} \times \mathbb{R}^N \times \boldsymbol{\Theta} \mapsto\boldsymbol{\mathcal{M}}$ is a Lipschitz function. The subscript notation $t|t-1$ is used to emphasize the fact that $\boldsymbol{\mu}_{t|t-1}$ is an approximation of the dynamic location process at time $t$ given the past, that is equivalent to say that $\boldsymbol{\mu}_{t|t-1}$ is $\mathcal{F}_{t-1}$-measurable. Therefore, based on past observations and a starting value $\boldsymbol{\mu}_{1|0} \in \boldsymbol{\mathcal{M}}$, one can approximate the unobserved path of $\boldsymbol{\mu}_t$ in (ref) by mimicking the recursion in (ref). It is typically assumed that a parameter value $\boldsymbol{\theta}_0$ exists, at which the true location can be recovered, i.e. $\boldsymbol{\mu}_{t|t-1}(\boldsymbol{\theta}_0) = \boldsymbol{\mu}_t$ (assumption (ref) of correct specification).
In this paper, we approximate the temporal changes of the dynamic location by relying on the score-driven framework of CrealKoopLucas2013 and Harvey2013. Specifically, we assume that, conditional on the past, the distribution of $\boldsymbol{y}_t$ is Student's t, with $\nu>0$ degrees of freedom and conditional location equal to $\boldsymbol{\mu}_{t|t-1}$, i.e.
and specify the SRE in (ref) as follows,
where $\boldsymbol{\omega}$ is a $\mathbb{R}^N$ vector of unconditional means, $\boldsymbol{\Phi}$ and $\boldsymbol{K}$ are $\mathbb{R}^{N \times N}$ matrices of coefficients and the driving force $\boldsymbol{u}_t$ is proportional to the score of conditional density in (ref). Indeed, the conditional score with respect to the time varying location filter is
where
with $w_t = 1 + (\boldsymbol{y}_t - \boldsymbol{\mu}_{t|t-1})^\top \boldsymbol{\Omega}^{-1} (\boldsymbol{y}_t - \boldsymbol{\mu}_{t|t-1})/\nu$, is a martingale difference sequence, i.e. $\mathbb{E}_{t-1}[\boldsymbol{u}_t]=\boldsymbol{0}_N$, under correct specification, where the shorthand notation $\mathbb{E}_{t-1}[X]$ is used for the conditional expectation $\mathbb{E}[X | \mathcal{F}_{t-1}]$. The score as the driving force in an updating equation for a time varying parameter is the key feature of score-driven models. The rationale behind the recursion (ref) is very intuitive. Analogously to the Gauss-Newton algorithm, it improves the model fit by pointing in the direction of the greatest increase of the likelihood. Optimality of score driven updates in observation-driven models is discussed by {Blasques_Koopman_Lucas2015.
In the context of location estimation under the Student's t assumption, a further relevant motivation for the score-driven methodology lies in the robustness of the implied filters. Indeed, the positive scaling factors $w_t$ in equation (ref) are scalar weights that involve the Mahalanobis distance. They possess the role of re-weighting the large deviation from the mean incorporated in the innovation error
Robustness comes precisely from winsorizing the innovation error $\boldsymbol{v}_t$. Note that when $\nu\rightarrow \infty$, $\boldsymbol{u}_t$ converges to $\boldsymbol{v}_t$ and equations ((ref)) and ((ref)) coincide with the steady state innovation form of a linear Gaussian state-space model.
A formal proof of the robustness of the method is in the following Lemma, which provides sufficient conditions for a filter to be robust, in line with Calvet_Czellar_Ronchetti2015. We first enounce the correct specification assumption.
The moment structure reveals important features of the driving force $\boldsymbol{u}_t$, that turns out to be an an IID sequence with zero mean vector and $(\operatorname{vec})$-variance covariance matrix,
Let us combine equations (ref) and (ref) and write the filter explicitly, as follows,
By starting at some initial value, $\boldsymbol{\mu}_{1|0}\in \boldsymbol{\mathcal{M}}$, and using equation (ref) for $t=1, \dots, T$, with $T \in \mathbb{N}$, one can recover a unique filtered path $\{ \hat{\boldsymbol{\mu}}_{t|t-1} \}_{t \in \mathbb{N}}$ for every $\boldsymbol{\theta} \in \boldsymbol{\Theta}$. A desirable property is that the values used to initialise the process are asymptotically negligible, in the sense that as the time $t$ increases, the impact of the chosen $\boldsymbol{\mu}_{1|0}$ eventually vanishes and the process will converge to a unique stationary and ergodic sequence. This stability property of the filtered sequence $\{ \hat{\boldsymbol{\mu}}_{t|t-1}\}_{t \in \mathbb{N}}$ is known as invertibility, see Straumann_Mikosh2006 and Blasques_Gorgi_etal2018. Existence of the unique stationary and ergodic solution to the SRE (ref) is established by Lemma (ref). Invertibility of the filter is proved in Lemma (ref).
The stability condition $\varrho(\boldsymbol{\Phi}) < 1$ is a well-known condition in the theory of linear systems, see Hannan1970, Hannan_Deistler1987 or Lutkepohl2007, which extends to the case of the present nonlinear model.
With the next Lemma, the relevant conditions under which the SRE in (ref) is contractive on average are given so that the convergence of the filtered sequence $\{ \hat{\boldsymbol{\mu}}_{t|t-1} \}_{t \in \mathbb{N}}$ to a unique $\mathcal{F}_{t-1}$-measurable stationary and ergodic solution $\{ \tilde{\boldsymbol{\mu}}_{t|t-1} \}_{t \in \mathbb{Z}}$, irrespective of the initialization $\boldsymbol{\mu}_{1|0}$, is obtained as a corollary of Theorem 3.1 of Bougerol1993 or, equivalently, of Theorem 2.8 of Straumann_Mikosh2006. Moreover, as a consequence of Lemma (ref), both $\{ \hat{\boldsymbol{\mu}}_{t|t-1} \}_{t \in \mathbb{N}}$ and $\{\tilde{\boldsymbol{\mu}}_{t|t-1} \}_{t \in \mathbb{Z}}$ have bounded moments.
The contraction condition in equation (ref) imposes restrictions on the parameter space $\boldsymbol{\Theta}$ that cannot be checked directly. Also, the expectation in the same equation cannot be verified in practice, since it depends on the unconditional, unknown, distribution of $\boldsymbol{y}_t$, see also the discussion in Blasques_Gorgi_etal2018. Thus, one can rely on sufficient conditions which are typically more restrictive than (ref) and that we discuss in the following, similarly to Linton2020. Specifically, the contraction condition in (ref) is satisfied if
Motivated by Example 3.8 of Straumann_Mikosh2006, we rewrite $\boldsymbol{X}_1$ at $\boldsymbol{\theta}_0$, so that equation (ref) becomes
Since $ \boldsymbol{\epsilon}_1 \sim \boldsymbol{t}_{\nu_0} (\boldsymbol{0}_N, \boldsymbol{I}_N)$, based on Monte Carlo simulations, Figure (ref) displays a region for a bivariate model ($N=2$) that satisfies the condition (ref) on a grid of values $(\|\boldsymbol{\Phi}_0\|, \|\boldsymbol{K}_0\|)\in (0,1)^2$, with $\nu_0 = 7$ and $\boldsymbol{\Omega}_0=\boldsymbol{I}_2$.
As expected, the restrictions that need to be imposed on $\| \boldsymbol{\Phi}_0 \|$ and $\| \boldsymbol{K}_0 \|$ are always stronger than those required for strict stationarity and ergodicity, see Lemma (ref). Neverthless, the region depicted in Figure (ref) shows that a subset $\boldsymbol{\Theta^*}$ of the parameter space $\boldsymbol{\Theta}$ exists, with $\|\boldsymbol{\Phi}\| < 1 $ and $\|\boldsymbol{K}\|$ sufficiently small such that (ref) is satisfied $\forall \boldsymbol{\theta} \in \boldsymbol \Theta^* \subset \boldsymbol{\Theta}$, producing a non degenerate invertibility region.
In alternative to the simulation-based method, one can restrict the estimation procedure to the empirical version of the invertibility constraint in (ref) as in Wintenberger2013 and Blasques_Gorgi_etal2018. The empirical counterpart of (ref) with $k=1$ is
for some $\delta>0$ arbitrarily small.
To conclude, we note that the process $\{ \boldsymbol{y}_t \}_{t\in\mathbb{Z}}$ inherits some properties from those of the filter evaluated at the true parameter value. As a consequence of Lemma (ref) and Lemma (ref), we obtain the following result.
Finally, the multi-step forecasts can be straightforwardly obtained as
Let ${\ell}_t(\boldsymbol{\theta})$ denote the conditional $\log$-likelihood function for a single observation, obtained by taking the logarithm of (ref) considered as a function of the parameter $\boldsymbol{\theta} = (\boldsymbol{\xi}^\top, \boldsymbol{\psi}^\top)^\top \in \boldsymbol{\Theta} \subset \mathbb{R}^p$, $\boldsymbol{\xi} = (\nu, (\operatorname{vech}(\boldsymbol{\Omega}))^\top, \boldsymbol{\omega}^\top)^\top \in \mathbb{R}^{s}$, with $s = 1 + \frac{1}{2}N(N+1) + N$ and $\boldsymbol{\psi} = ( (\operatorname{vec}\boldsymbol{\Phi})^\top, (\operatorname{vec}\boldsymbol{K})^\top)^\top\in \mathbb{R}^{d}$, with $d = (N \times N) + (N \times N)$ and hence, $p = s + d$.
Lemma (ref), ensures that any choices of the initial condition $\boldsymbol{\mu}_{1|0} \in \boldsymbol{\mathcal{M}}$ used for starting the filtering process are asymptotically equivalent, such that, once an initial value has been fixed, it is possible to obtain an approximated version of the conditional $\log$-likelihood, $\hat{\ell}_t(\boldsymbol{\theta})$, by replacing $\boldsymbol{\mu}_{t|t-1}$ in ${\ell}_t(\boldsymbol{\theta})$ by the filtered dynamic location $\hat{\boldsymbol{\mu}}_{t|t-1}$. Thus, for the whole sample, we obtain $\hat{\ell}_T(\boldsymbol{\theta}) = \sum_{t=1}^{T} \hat{\ell}_t(\boldsymbol{\theta})$ and the MLE of $\boldsymbol{\theta}$ is
We now discuss strong consistency and asymptotic normality of the MLE. The following assumptions are standard in the likelihood theory of non linear observation driven models.
Assumption $4.1.1$ can be replaced by the conditions of Lemma (ref). Assumption $4.1.2$ ensures that the filtered sequence $\{\hat{\mu}_{t|t-1}\}_{t\in\mathbb{N}}$ converges to a stationary ergodic limit sequence, irrespective of the initial conditions. Assumptions $4.1.3$ and $4.1.4$ ensure the existence of the MLE and the validity of first order asymptotics. Assumption $4.1.5$ guarantees the existence of the information matrix.
By Theorem (ref), $\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}_0)$ can be consistently estimated by
As the dynamic location and its derivatives are nonlinear functions of the parameter $\boldsymbol{\theta}$, the general formula for the second derivatives in (ref) has the form below
To avoid the recursive evaluation of the second derivatives of the dynamic location vector, a simpler consistent estimator can be obtained based on the analytical form of the conditional information matrix $\boldsymbol{\mathcal{I}}_{t}(\boldsymbol{\theta})$, as in Fiorentini_Sentana_Calzolari2003, defined as
Indeed, by the law of iterated expectations, one has
Given the assumption of correct specification, the score vector evaluated at the true parameter vector $\boldsymbol{\theta}_0$ forms a martingale difference sequence, so that, under the assumptions of Theorem (ref), asymptotic results for martingale difference sequences can be applied. In addition, the dynamic location (and its derivatives) are $\mathcal{F}_{t-1}$-measurable functions and therefore, after taking the conditional expectation, the last term in the right-hand-side of equation (ref) will cancel out.
It follows that, by Theorem (ref), $\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}_0)$ can be consistently estimated by
where $\widehat{\boldsymbol{\mathcal{I}}}_t( \widehat{\boldsymbol{\theta}}_T) $ is the conditional information matrix in (ref) evaluated at the filtered dynamic location $\hat{\boldsymbol{\mu}}_{t|t-1}$ and at the MLE $\widehat{\boldsymbol{\theta}}_T$. The analytical form of $\boldsymbol{\mathcal{I}}_{t}(\boldsymbol{\theta})$ is derived in section S2.3.
ML estimation and inference are carried out by means of Fisher's scoring method. A strongly reliable algorithm based on analytical formulae for the score vector and the Hessian matrix (reported in Appendix S2) is developed, which can be directly implemented in any statistical package through the following steps:
The analytical expressions for the score vector and the conditional information matrix used in step $2$ are in Section S2.
To initialise the estimation procedure, we follow the approach suggested in Fiorentini_Sentana_Calzolari2003. First, a consistent estimator of the restricted version of the parameter vector $\tilde{\boldsymbol{\theta}}_T$ is obtained by the Gaussian quasi-ML procedure in Bollerslev_Wooldridge1992. Second, a consistent method of moments is adopted for the degrees of freedom $\nu$, by making use of the empirical coefficient of excess kurtosis $\tilde{\kappa}$ on the standardized residuals and of the relation $\tilde{\nu} = (4\tilde{\kappa} + 6)/\tilde{\kappa}$. Convergence is fast in that usually few iterations of that procedure are needed, which makes scoring methods particularly appealing for estimation purposes.
In section S1, we report the details of a Monte Carlo study aimed to assess the finite sample properties of the MLE based on the Fisher's scoring method detailed in the above section. In summary, our approach performs well in terms of bias and root mean square errors for a wide range of time series, from the most severe heavy-tailed case (i.e., $\nu$ very small) to the Gaussian case (i.e. for $\nu\rightarrow \infty$), thus covering also the case of potential misspecification. In addition, we note that it delivers satisfactory results even when the number of iterations of the algorithm is limited to ten rounds.
In order to demonstrate a potential use of the robust score-driven filter, we show an innovative application to the estimation of consumer prices from homescan data. This field of application is gaining interest, due to the growing availability of high frequency and high detail purchase data collected through scanner technologies at the retail point (retail scan) or household level (homescan). The latter of type of data allows one to obtain cost-of-living measures for vulnerable sub-groups of the population, and to explore the distributional effects of fiscal measures. While being a valuable source for detailed price information, post-purchase homescan price data are affected by a measurement noise that can be potentially large in small samples, and the application of filtering techniques may help to mitigate such noise and control for outliers.
Scanner data are collected either at the retail level, e.g. supermarket data, or from households in consumer panels, i.e. homescan data. Retail scanner data are widely used to estimate prices, both for continuity with the traditional price survey methodology, and because they are expected to suffer less from the substitution (unit value) bias (Silver2001). This bias is due to the fact that scanner data are based on actual transactions, i.e. prices are only observed after the consumer purchases the good. This implies that the observed price embodies a quality choice component, as consumers confronted with a price increase may opt for a cheaper option (or a cheaper retailer) and information on non-purchased items is missing. The bias can be particularly important for aggregated goods, such as those goods commonly represented by category-level prices like food and drinks. Thus, a wide body of research has been devoted to improve sampling strategies and the choice of weights in aggregation. A well-documented problem is the change in the composition of the consumption basket over time, an issue that can be exacerbated by high-frequency data Feenstra2003b. For example, stockpiling of goods during promotion periods generate bias in price indices, as the purchased quantities are not independent over subsequent time periods Ivancic2011,Melser2018.
Although supermarket-level scanner data allow to mitigate the problem, as one expects a wide range of products to be purchased across the population of customers within a given time period, the use of homescan data to estimate prices and price indices has potentially major advantages. These advantages lie in the possibility to exploit household-level heterogeneity. Most importantly, it becomes feasible to estimate prices faced by particular population sub-groups whose consumption basket differs from the average one, as elderly households or low-income groups Kaplan2017,Broda2009. However, the unit value issue is heavier with homescan data, as individual households buy a small range of products. Thus, variable shopping frequencies and zero purchases make it necessary to rely on very large samples of households to control the bias. The problem becomes even more conspicuous for prices at the regional level, for products that are not frequently purchased and for products whose demand is highly seasonal.
Robust filtering techniques may constitute a powerful solution to the above mentioned problems, and may perform well even with relatively small samples of household as the one used in our application.
To illustrate the potential contribution of the proposed method, we exploit a data-set that has been recently used to evaluate the effects of a tax on sugar-sweetened beveraged introduced in France in 2012 Capacci2019. Our data consists of weekly scanner price data for food and non-alcoholic drinks. The data were collected in a single region, within the Italian GfK homescan consumer panel, based on purchase information on 318 households surveyed in the Piedmont region, over the period between January 2011 and December 2012. The regional scope and the relatively small sample provide an ideal setting to test the applicability and effectiveness of the multivariate filtering approach.
The data for our application consist of three time series of weekly unit values for food items, non-alcoholic drinks and Coca-Cola purchased by a sample of 318 households residing in the Piedmont region, Italy, over the period 2011-2012, and collected within the GfK Europanel homescan survey. The data-set provides information on weekly expenditures and purchased quantities for each of the three aggregated items, and unit values are obtained as expenditure-quantity ratios.
Average unit values are shown in Table (ref). Food and non-alcoholic drinks are composite aggregates, hence they are potentially subject to fluctuations in response to changes in the consumer basket even when prices are stable. Instead, Coca-Cola is a relatively homogeneous good, with little variability due to different packaging sizes.
We fit the multivariate score-driven model developed in the paper to the considered vector of time series. ML estimation produces the following multivariate dynamic system of time varying locations for Drinks (D), Food (F) and Coca-Cola (C),
where the values in parenthesis are the standard errors and with \[ \hat{\nu}= 6.921 \,\,\,{(0.229)} , \hat{\boldsymbol{\Omega}}=
\times 10^{-3}. \]
The estimated degrees of freedom are approximately $7$. We remark that the assumption of a (conditional) multivariate Student's t distribution implies that all the univariate marginal distributions are tail equivalent, see Resnick2004. This requires the implicit underlying assumption that the level of heavy-tailedness across the observed time series vector is fairly homogeneous. To investigate this issue, and for the sake of comparisons, we have carried out a univariate analysis, as in HarveyLuati2014, from which it resulted that the estimated degrees of freedom were very low for Coca-Cola (about $4$) and medium size (smaller than $30$) for the other two series, as expected. Hence, the multivariate score-driven model developed in the paper reveals to be a good compromise between a multivariate non-robust filter, based on a linear Gaussian model, and a robust univariate estimator. Indeed, a multivariate Portmanteau test on the residuals obtained from the three univariate models is carried out to test the null hypothesis $H_0 : \boldsymbol{R}_1 = \dots = \boldsymbol{R}_m = \boldsymbol{0}$, where $\boldsymbol{R}_i$ is the sample cross-correlation matrix for some $i \in \{1, \dots, m\}$ against the alternative $H_1 : \boldsymbol{R}_i \neq \boldsymbol{0}$. The results of Table (ref) indicate rejection of the null hypothesis of absence of of serial dependence in the trivariate series at the $5\%$ significance level.
We also remark that the estimated degrees of freedom close to $7$ rule out the hypothesis that the data come from a linear Gaussian state-space model, in which case the estimated degrees of freedom would be definitely higher. Nevertheless, we have fitted a misspecified linear Gaussian state-space model estimated with the Kalman filter and, as expected, along with a higher sensitivity to extreme values, in particular in the last period of the Coca-Cola series, likelihood and information criteria are in favour of the multivariate model based on the conditional Student's $t$ distribution.
The matrix of the estimated autoregressive coefficients $\hat{\boldsymbol{\Phi}}$ measures the dependence across the filtered dynamic locations $\hat{\boldsymbol{\mu}}_{t|t-1}$, while the estimated scale matrix $\hat{\boldsymbol{\Omega}}$ measures the concurrent relationship between the three series under investigation, i.e. drink, food and Coca-Cola prices. For these matrices, we report the estimates of the coefficients and, in parenthesis, the relative standard errors. The diagonal elements of $\hat{\boldsymbol{\Phi}}$ show that each variable of interest is highly persistent. In order to explore the relation among the series, we implement an impulse response analysis. Figure (ref) shows the estimated impulse response functions.
The nonlinear impulse are computed by using the local projections approach of Jorda2005, and the confidence bands are obtained by using the Newey-West corrected standard-errors, see Newey1987. What emerges is a negative relation between drink and food prices: a unit shock in drink prices will produce a negative shock in food prices. This may adjustments in purchasing decisions by the households aimed at mitigating the rising cost of their shopping basket. This would be evidence that univariate signals are likely to suffer from the unit value bias. Similarly, a non trivial negative relation exists between food and Coca-Cola prices. A unit shock on food prices yields a concurrent negative impact on Coca-Cola prices, which is also noted from the analysis of the cross-correlations. As one might expect, a positive correlation exists between Coca-Cola prices and drink prices, as the former product belongs to the latter category. Instead, unit shocks on food prices seem to have negligible correlation (if any) on drink prices.
Figure (ref) shows the original unit value time series and the corresponding signals extracted through the multivariate score-driven filter. Noise and outliers, as well as some irregular periodic pattern, are clearly visible in the drinks and food series. On the other hand, the Coca-Cola series is relatively regular, with the exception of few peaks, including a couple of large outliers in the second year. Given the homogeneous nature of the good, it is reasonable to believe that those extreme values are the results of measurement error.
The estimates illustrate an effective noise reduction and return patterns that are smoother and more consistent with a regular price time series. As one would expect, the Coca-Cola DCS-t series is very flat, and suggests a relatively stable price over the two-years time window, with no outliers.
Figure (ref) shows the monthly natural logarithm differences of the raw homescan prices (HSP) and the estimated signals, together with changes in the official Regional CPIs (R-CPI) for food and non-alcoholic drinks, whereas no CPI to the brand detail is produced. The R-CPIs are provided by the National Statistical Institute (ISTAT). They have a monthly frequency and are built with a traditional survey-based approach on retailers. The comparison between the score-driven filtered values and the R-CPIs is purely indicative, as the unit values from the homescan data are weekly, whereas the official CPIs are monthly. This frequency difference may lead to biased comparisons Diewert2016. Nevertheless, the graphs confirm that the score-driven signals are effective in reducing the noise in the data. This is especially true for the food series, whose CPIs are more volatile compared to drinks. The correlation between the raw homescan log-differenced unit value and the log-differenced food CPI is 0.05, against 0.44 when the filtered time series is considered. For the non-alcoholic drinks price series the gain is less conspicuous, as prices evolve very regularly over the time window. Still, an inexistent correlation between the HSP and the R-CPI (-0.02) turns into a positive one (+0.11) when considering the score-driven estimates and the R-CPI.
In essence, the empirical evidence suggests that a robust multivariate approach to model-based signal extraction produce meaningful price series from homescan data, especially when noise and outliers in the original data are relevant. We find the approach to perform reasonably well even with a low number of sampled households (318) and price time series (3), and with a relatively short time window (104 weeks). Future research might shed further light on the implications of dealing with a larger number of price series and longer time series.
We developed a nonlinear and multivariate dynamic location filter which enables the extraction of reliable signals from vector processes affected by outliers and possibly non-Gaussian errors. Its peculiarity lies in the specification of a score-robust updating equation for the time-varying conditional location vector. Compared to the existing literature on observation driven models for time varying parameters, the model has two innovative features: (a) it extends the univariate first-order dynamic conditional location score by HarveyLuati2014 to the multivariate setting; and (b) it extends the dynamic model for time varying volatilities and correlations by CrealKoopLucas2011 to the location case.
We derived the stochastic properties of the filter and, under correct specification, of the data generating process: bounded moments, stationarity, ergodicity, and filter invertibility. Parameters are estimated by ML and we provided closed formulae for the score vector and the Hessian matrix, which can be directly used for a scoring procedure. Consistency and asymptotic normality have been proved and a Monte-Carlo study showed good and reliable finite sample properties. In the case when the degrees of freedom tend to infinity, or, in practice, their estimate is of the order of hundreds, our specification converges to a linear and Gaussian model.
The empirical application showed that robust filtering may lead to satisfactory estimates of price signals from homescan data, in the case when the multivariate dimension is low. We contribute to research in this area with two promising results. First, we show that robust modeling allowing for heavy tails is more effective in dealing with noisy series affected by outliers or extreme observations. Second, the multivariate extension of the DCS-t model has shown more appropriate than the robust univariate filtering approach in the case of scanner price data, as price time series are expected to have a good degree of correlation. This proves to be valuable information to reduce the noise across the modelled price time series.
Additional supporting information may be found in the online appendix for this article at the publisher's website.