A Comment on ”The Art of Data Augmentation”

Gabriel Huerta,

Department of Probability and Statistics, Centro de Investigación en Matemáticas,

Apartado Postal 402, Guanajuato, Gto. 36000, México

&

Wenxin Jiang and Martin A. Tanner^{1}

Department of Statistics, Northwestern University, 2006 Sheridan Rd., Evanston IL 60208-4070, U.S.A.

1M. Tanner was supported by NIH grant CA35464.

This insightful paper by van Dyk and Meng (VM) makes the important point that advances in data augmentation algorithms offer a wide variety of tools for statistical inference. Time series methods are no exception and mixture modeling within this context may help to improve forecasting and to detect changes in structure across time.

The time series modeling approach that we adopt is based on the idea of mixing models through the neural network paradigm known as Hierarchical Mixtures-of-Experts (HME)

- see Jordan and Jacobs (1994). The HME approach easily allows for model comparison and permits one to represent the mixture weights as a function of time or other covariables. With the additional hierarchy, it is possible to localize the comparisons to specific regions or regimes. Furthermore, the defining elements of the mixture do not have to be restricted to a particular class of models permitting very general comparisons. In this comment parameters are estimated via maximum likelihood using the EM-algorithm- extensions to a full Bayesian approach using MCMC may follow one or more of the many lines outlined by VM. We see this comment as a call to the Chagalls of this world to use their artist abilities to develop quick mixing stochastic algorithms for this important, yet complex class of HME models.

Let {y_{t}}_{0}^{n} be a time series of endogenous or response variables, and {x_{t}}_{0}^{n} be a time series of exogenous variables or covariates. Suppose the main interest is to draw inference on {y_{t}}_{0}^{n} conditional on {x_{t}}_{0}^{n}. Let the conditional probability density function (pdf) of y_{t} be f_{t}(y_{t}|F_{t−1}, X; θ), where θ is a parameter vector; X is the σ-field generated by {x_{t}}_{0}^{n}, representing the external information; and for each t, F_{t−1} is the σ-field generated by {y_{s}}_{0}^{t−1} representing “the previous history” at time t − 1. Typically, the conditional pdf f_{t} is assumed to depend on X through x_{t} only. In HME, the pdf f_{t} of the response variable is assumed to be a conditional mixture of the pdfs from simpler, well established models. In a time series context, this mixture can be represented by the finite sum

f_{t}(y_{t}|F_{t−1}, X; θ) = ∑_{J} g_{t}(J|F_{t−1}, X; γ)π_{t}(y_{t}|F_{t−1}, X , J; η),  (1)

where the functions g_{t}(·|·, ·; γ) are the mixture weights; π_{t}(·|·, ·, J; η) are the conditional pdfs from simpler models each defined by a label J; and γ and η are vectors of sub-parameters from θ.

The simpler models in HME are often referred to as the “experts”. In a time series context, one “expert” could be an AR(1) model, another “expert” could be a GARCH(1,1) model or an EGARCH(1,1) model. For example, in a situation where it is not clear whether to use a stochastic or a deterministic trend, one expert could be a trend-stationary process, another a difference-stationary process. A somewhat simpler situation occurs when all the experts propose a model of the same type, e.g. linear autoregressive, but perhaps with different values for the coefficients or for the model order.

Furthermore, the HME models considered have an additional layer designed with the purpose of local time series modeling. The HME partitions the covariate space, which could include time, into O overlapping regions called “overlays”. In each overlay, M models are to compete with each other, in the hope that the model most suitable to the specific region is favored by a high weight (see Figure 1). By having multiple overlays, the hierarchical mixture model allows for modeling multiple switching across regions.

Therefore, the expert index J can be expressed as J = (o, m), where the overlay index o takes a value from {1, …, O}, and the model-type index m from {1, …, M}. We allow the same type of model m to assume different versions or more specifically different parameter values, at each possible overlay.

The mixing weights are often referred to as “gating functions”. They can depend on the previous history, exogenous information (see McCulloch and Tsay, 1993), or can exclusively depend on t. The gating functions may have the form

g_{t}(o, m|F_{t−1}, X; γ) = \{\frac{e^{v_{o}+u_{o}^{T} w_{t}}}{∑_{s=1}^{O} e^{v_{s}+u_{s}^{T} w_{t}}} \} \{\frac{e^{v_{m|o}+u_{m|o}^{T}w_{t}}}{∑_{l=1}^{M} e^{v_{l|o}+u_{l|o}^{T}w_{t}}} \} ,  (2)

where the v’s and u’s are parameter components of γ; and w_{t} is an “input” at time t, which is measurable with respect to the σ-field induced by F_{t−1} ∪ X. For example, the input w_{t} could be the covariate x_{t}, the “two-lag” history (y_{t−1}, y_{t−2})^{T}, or exclusively depend on time t.

In the context where one is interested in how the weighting for individual models is assigned across different time periods, w_{t} can be taken as (t/n). Therefore, one can adopt the following parametric form for the gating functions:

g_{t}(o, m|F_{t−1}, X; γ) = g_{om}(t; γ) ≡ \{\frac{e^{v_{o}+u_{o}(t/n)}}{∑_{s=1}^{O} e^{v_{s}+u_{s}(t/n)}} \} \{\frac{e^{v_{m|o}+u_{m|o}(t/n)}}{∑_{l=1}^{M} e^{v_{l|o}+u_{l|o}(t/n)}} \} . (3)

Here γ includes all the following components: v_{1}, u_{1}, …, v_{O−1}, u_{O−1}, …, v_{1|1}, u_{1|1}, …, v_{M−1|1}, u_{M−1|1}, …, v_{M−1|O}, u_{M−1|O}. For identifiability, we set v_{O} = u_{O} = v_{M|o} = u_{M|o} = 0 for all o = 1, …, O. The free vector of parameters γ in the gating functions automatically determines the location and the “softness” of the splitting of the regions.

Note that this framework defines the two-layer HME architecture of Jordan and Jacobs (1994), where the first layer of gating functions hypothesizes O overlays on the entire time axis, and the second layer of gating functions defines weights for each of the M model types within each overlay. When the input space for the gating functions is time, the hierarchical mixture model can identify the region over which a model or a set of models is (are) dominant in a data-adaptive manner. Thus, the present approach allows for modeling multiple regime switching. Further details of this approach, as well as related asymptotic theory are presented in Huerta, Jiang and Tanner (2000).

Figure 1: A Graphical Representation of a Two-Layer HME.

Inference on the parameter θ can be based on the log-likelihood function, conditional on y_{0}, X and “averaged” in time, which is

L_{n}(·) = n^{−1} ∑_{t=1}^{n} log f_{t}(y_{t}|F_{t−1}, X; ·).  (4)

We denote the maximum likelihood estimate (MLE) of θ as θ̂ = arg max L_{n}(·). To obtain the MLE, the EM algorithm starts with an initial estimate of the parameters θ^{0}. Then a sequence {θ^{i}} is obtained by iterating between the following two steps:

For i = 0, 1, 2, …,

\underline{E-step:} Construct

Q^{i}(θ) = ∑_{t=1}^{n} ∑_{o,m} h_{om}(t; θ^{i}) log{π_{t}(y_{t}|F_{t−1}, X , o, m; η)g_{t}(o, m|F_{t−1}, X; γ)}, (5)

where θ = (γ, η), θ^{i} = (γ^{i}, η^{i}), h_{om}(t; θ^{i}) = h_{om}(t; θ)|_{θ=θ^{i}}, and

h_{om}(t; θ) = \frac{g_{om}(t; γ)π_{t}(y_{t}|F_{t−1}, X , o, m; η)}{∑_{s=1}^{O} ∑_{l=1}^{M} g_{sl}(t; γ)π_{t}(y_{t}|F_{t−1}, X , s, l; η)}  (6)

is the “posterior probability” of choosing the expert (o, m) at time t.

\underline{M-step:} Find θ^{i+1} = arg max_{θ} Q^{i}(θ).

Inference is greatly facilitated by the introduction of augmented data, resulting in the fact that the objective function Q^{i} has the form of a double sum of logarithms, instead of a “sum log sum” typical for the log likelihood function L_{n}. For this reason, the maximization of the objective function can be decomposed into a number of smaller maximization problems which involve fewer parameters and usually define “known” maximizations of widely used models. For example, suppose the expert pdf has the form

π_{t}(y_{t}|F_{t−1}, X , o, m; η) = p_{t}(y_{t}|F_{t−1}, X , m; η_{om}),  (7)

where η is decomposed into a collection of sub-parameter η_{om}, each of which only appears in the pdf of one expert. The parameter η_{om} carries an index o in addition to m to allow one type of model to take different versions (parameters) in different overlays. In such a situation, in the M step, the maximization over the η_{om}’s and γ can be performed separately. For example, for each o, m, i,

η_{om}^{i+1} = arg max_{η_{om}} ∑_{t=1}^{n} h_{om}(t; θ^{i}) log p_{t}(y_{t}|F_{t−1}, X , m; η_{om}),  (8)

which become the “standard” (albeit weighted by the h’s) maximum likelihood problem for model type-m.

When the MLE θ̂ is obtained, we are interested in evaluating the relative weighting of each of the M model types at time t. Two estimates are of interest in this regard. One is an empirical Bayes estimate of the posterior probability / weight of model type-m at time t. This is the conditional probability regarding the history up to time t and defined by:

P̂_{t}(m|y_{t}, F_{t−1}, X) ≡ ĥ_{m}(t) ≡ ∑_{o=1}^{O} h_{om}(t; θ̂),  (9)

where θ̂ is the MLE. Another approach for weighting is to consider an empirical Bayes-type estimate of the unconditional probability / weight of model m at time t (unconditional on the history of the endogenous process {y_{t}}):

P̂_{t}(m|F_{t−1}, X) ≡ ĝ_{m}(t) ≡ ∑_{o=1}^{O} g_{om}(t; θ̂).  (10)

As we shall see in an example, (9) can vary point-wise over time due to the conditioning on the specific history of the observations. The second weighting scheme (10) is smoother when describing a regional change of preference for model m. The term ĥ_{m}(t) is an estimated “posterior” probability of model m, and ĝ_{m}(t) is the corresponding estimate of the “prior” probability in the sense that the prior probability is not conditional y_{t} and the posterior probability is conditional on y_{t}. These estimates are not formal Bayesian priors and posterior probabilities for model m, since we have not assigned any prior p(θ) to the parameters θ, but instead are estimating the conditional or “unconditional” mixing weights at the MLE.

Figure 2: Exchange Rate Data. Top: 500 observations of spot rates between the Germanmark and US-dollar beginning in October 1986. Bottom: Logarithm of the first difference of the spot rates between the German-mark and the US-dollar.

We consider a financial time series to illustrate this hierarchical mixture defined with GARCH and EGARCH models - see Bollerslev (1986) and Nelson (1991). The series consist of 500 daily observations of exchange rates between the German-mark and the US-dollar starting from October 9, 1986. In fact, Figure 2 presents a time plot of the 500 daily spot rates characterized by a non-stationary random walk behavior.

Also in the same figure, we present the logarithm of the first difference of spot rates as a function of time, a transformation which is widely used to induce covariance-stationarity and propose parametric models. Assuming that Y_{t} represents the log of the first difference in spot rates, as discussed in Andersen and Bollerslev (1998), 3 candidate models are used to obtain inferences on volatilities or innovation variances at each time t:

• An AR(1) model simply defined by:

Y_{t} = ϕ_{0} + ϕ_{1}Y_{t−1} + ϵ_{t},

where ϵ_{t} ∼ N(0, σ^{2}).

• An AR(1)-GARCH(1,1) model represented by

Y_{t} = ϕ_{0} + ϕ_{1}Y_{t−1} + ϵ_{t},

but now ϵ_{t} ∼ N(0, σ_{t}^{2}). That is, the innovation variance can change in time and according to the evolution equation

σ_{t}^{2} = θ_{0} + θ_{1}ϵ_{t−1}^{2} + θ_{2}σ_{t−1}^{2},

with non-negative parameters θ_{0}, θ_{1} and θ_{2}.

• Finally, an AR(1)-EGARCH(1,1) is considered, where again

Y_{t} = ϕ_{0} + ϕ_{1}Y_{t−1} + ϵ_{t}

but now the evolution in variance is defined in terms of the natural logarithm and the standardized innovations z_{t} = ϵ_{t}/σ_{t} through the expression

log(σ_{t}^{2}) = β_{0} + β_{1}z_{t−1} + β_{2}(|z_{t−1}| − \sqrt{2/π}) + β_{3} log(σ_{t−1}^{2}),

with no restrictions on the parameters β_{0}, β_{1}, β_{2} and β_{3}.

Within an HME framework taking M = 3: π_{t}|_{o,m=1} denotes the pdf of Y_{t} given the history based on an AR(1); π_{t}|_{o,m=2} the same pdf with an AR(1)-GARCH(1,1); and π_{t}|_{o,m=3} denotes the pdf with an AR(1)-EGARCH(1,1). As before, the index o is added to model parameters and our initial exploration is based on a value of O = 2, i.e. allowing for 2 overlays. We ran the EM-algorithm with 20 different starting points. Parameters for the pdfs π_{t}|_{o,m} were initialized at the individual model MLE’s and initial parameters for the gating functions were generated from uniform distributions. Each EM was run for 500 iterations and solutions were ranked using the log-likelihood function L_{n}(·).

Figure 3 presents the estimates ĝ_{m}(t) for m = 1, 2, 3. In general, the model assigns a very low weight to the AR(1), with competing weights for AR(1)-GARCH(1,1) and AR(1)- EGARCH(1,1), for approximately the first 100 observations of the series. For the remaining segment of the time series, the preferred model is the AR(1)-EGARCH(1,1).

In Figure 4, we present the estimates of ĥ_{m}(t) for m = 1, 2, 3. Although the general “smoothed” pattern is similar to that exhibited by Figure 3, the model posterior probabilities have jumps of high probability for the three competing alternatives. This example reflects how h_{m}(t) can be highly impacted by single observations. The “ups and downs” in volatility experienced by exchange rate data across time lead to these model switches. Periods of almost constant variance can be well represented by an AR(1) model but when the data present periods of non-constant variance, the GARCH or EGARCH structure dominates producing the large jumps in the functions h_{m}(t).

Figure 5 presents a comparison of the estimated volatilities of the HME with those based on the individual models. The MLE, computed with a numerical optimizer, was Figure 3: Exchange Rate Data Example. Maximum likelihood estimates of g_{m}(t) for the 3 models considered: AR(1), AR(1)-GARCH(1,1) and AR(1)-EGARCH(1,1).

used to obtain the volatilities for AR(1), GARCH(1,1) and EGARCH(1,1). For the present hierarchical mixture model, the volatilities were estimated using the EM-solution, recognizing that for this mixture model, the variance of Y_{t} given the history is the expectation with respect to the mixing weights g_{m}(t) of individual model-variances plus the variance of the expectation functions for each defining model. We note that the volatility for the HME smoothes some of the high volatility peaks induced by other models but recognizes the overall pattern suggested by the AR(1)-GARCH(1,1) and the AR(1)-EGARCH(1,1) processes.

Figure 4: Exchange Rate Data Example. Maximum likelihood estimates of h_{m}(t) for the 3 models considered: AR(1), AR(1)-GARCH(1,1) and AR(1)-EGARCH(1,1).

Figure 5: Exchange Rate Data Example. Comparison of estimated volatilities for the HME and the individual models AR(1), AR(1)-GARCH(1,1) and AR(1)-EGARCH(1,1).

REFERENCES

Andersen, T. G. and Bollerslev, T. (1998). ARCH and GARCH models. Encyclopedia of Statistical Sciences. Update Vol. 2, 6-16. Edited by Kotz, S., Read, C. B. and Banks, D. L., Wiley, New York.

Bollerslev T. (1986). Generalized autoregressive conditional heteroskedasticity Journal of Econometrics 31, 307-327.

Huerta, G., Jiang, W., and Tanner, M. A. (2000). Time series modeling via hierarchical mixtures. Technical Report, Northwestern University.

Jordan, M. I., and Jacobs, R. A. (1994). Hierarchical mixtures of experts and the EM algorithm. Neural Comp. 6, 181-214.

McCulloch, R. E. and Tsay, R. S. (1993). Bayesian inference and prediction for mean and variance shifts in autoregressive time series. Journal of the American Statistical Association 88, 968-978.

Nelson, D.B. (1991). Conditional heteroskedasticity in asset returns: a new approach. Econometrica 59, 347-370.