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.
91,979 characters · 16 sections · 40 citation commands
Analysis of Distributional Dynamics for Repeated Cross-Sectional and Intra-Period Observations
JEL classification codes: C13, C14, C22 \\ Key words and phrases: distributional dynamics, repeated panels, intra-period distributions, functional autoregression, time-varying density, density forecasting
\onehalfspacing
Many important issues in economics are related to the time evolution of state distributions of economic variables, which are defined to be either cross-sectional or intra-period distributions of the underlying economic variables. For example, the time series of cross-sectional or intra-period distributions of returns from traded assets contain important information on the dynamics of risk and return in financial markets. However, tracking the complete distributional dynamics is not an easy task because of the infinite dimensional nature of state distributions. The conventional approaches usually get around the problem by only considering some particular dimensions of the distributions such as the mean, the variance, or some specific quantiles. Information beyond these dimensions, however, would be necessarily ignored. In this paper, we develop a novel analytical framework that can be used to track and exploit the dynamics of entire state distributions.
We consider the functional autoregressive (FAR) model of time-varying distributions. Allowing the distributions to be infinite-dimensional, we adopt a nonparametric approach to modeling distributional dynamics. The parametric approach has been well studied, including the autoregressive conditional moment models by engle-82, bollerslev-86, engle-lilien-robins-87, harvey-siddique-99 and brooks-burke-heravi-persand-05, the autoregressive conditional density model by hansen-94, and the autoregressive quantile models by koenker-xiao-06 and xiao-koenker-09, among many others. Our approach, on the contrary, does not require any parametric assumptions on the underlying state distributions. Nevertheless, it allows us to investigate full temporal dynamics of state distributions, including dependence structures among moments and tail probabilities, within a single framework.
For our analysis, we represent state distributions by their density functions and regard them as random elements of a Hilbert space. Of course, it is also possible to use other functions such as distribution functions and quantile functions to represent state distributions.\footnote{boogaart-egozcue-pawlowsky-glahn-14 propose to use the so-called Bayes Hilbert space, which uses log-densities instead of densities themselves, to study distributional dynamics. Though their approach has some advantages, it does not allow for any of the methodologies developed in our paper to analyze distributional dynamics.} However, we choose density functions, since they provide clear and direct interpretations on the dynamics of moments, tail probabilities and other important features of state distributions. Indeed, in our framework, the cross-sectional moments and tail probabilities are defined simply as the inner products of state density functions with some deterministic functions in a Hilbert space. No other functions representing state distributions have this nice property. We fit our FAR model using the time series of state density functions, which are obtained from repeated cross-sectional and intra-period high frequency observations.\footnote{In all our applications, the number of cross-sectional or intra-period observations is relatively much bigger than the number of time series observations. Therefore, we ignore the statistical errors incurred by the estimation of state densities using repeated cross-sectional and intra-period observations.} Once fitted, we estimate the leading progressive/regressive features and conduct various impulse response analyses to reveal the way in which future state distributions respond to changes in past state distributions. In addition, we develop a variance decomposition method that relates a specific feature of a current state distribution to the moments of its past. Lastly, we propose a density predictor and find its asymptotic properties.
Recently, there has been a rapid development in the theory and applications of functional time series analysis. Many authors, including besse-cardot-stephenson-00, bosq-00, antoniaids-paparoditis-sapatinas-06, ferraty-vieu-06, mas-07, kargin-onatski-08, hormann-kokoszka-10, horvath-huskova-kokoszka-10, didericksen-kokoszka-zhang-12, kokoszka-reimherr-13, panaretos-tavakoli-13a, panaretos-tavakoli-13b, aue-norinho-hormann-15, aue-horvath-pellatt-17 and klepsch-kluppelberg-17, among others, studied various aspects of functional data. In many of these works, however, temporal dependencies in the functional dynamics are used in a relatively implicit way to make forecasts. We in this paper emphasize the investigation of the functional dynamics themselves in the distributional dynamics setting, and develop tools to interpret them in economically meaningful manners.
We apply our methodology to two representative state distributions: the intra-month distributions of the UK Pound/US Dollar (GBP/USD) exchange rate 15-minute log returns and the cross-sectional distributions of the New York Stock Exchange (NYSE) stocks monthly returns. In both cases, we find that observations/stocks with small returns play the most important role in interacting with the past distribution and in determining the future distribution. Outliers turn out to be unimportant in the distributional dynamics. We also find that various characteristics such as the moments and the tail probabilities of the current distribution respond differently to shocks to, and subsequent changes in, the relative frequency of observations with different levels of returns. The momentum effect appears in the dynamics of the first moment: The more observations with positive (negative) returns in the past distribution, the higher (lower) the first moment of the current distribution. On the other hand, sizes of the returns play a more important role in the dynamics of the second moment; the more observations with small (large) returns in the past distribution, the lower (higher) the second moment of the current distribution.
In contrast to the symmetry found in the behaviors of left and right tail probabilities in the foreign exchange rate return application, those in the stock return application appear to be quite asymmetric. In the former, a higher relative frequency of large returns in the past distribution, regardless of the signs of these returns, increases both the current left and right tail probabilities rather symmetrically. In the latter, however, the current left tail probability is mostly determined by the relative frequency of moderate-to-large negative returns in the past distribution, whereas the current right tail probability is positively affected both by the relative frequency of moderate-to-large negative returns, though to a lesser degree, and by the relative frequency of positive returns in the past distribution. The variance decomposition exercises indicate that in both cases, the variances of the current first moments are not explained much by the variances of the past integral moments, while the current second moments are more closely connected to the past second moments. These results lend support to the usual observation that in financial markets, the mean is not quite predictable while the volatility is much more persistent, and thus predictable. These moment dependence findings are also consistent with the rolling out-of-sample forecast results using our model.
The rest of the paper is organized as follows. In Section 2, we introduce the functional autoregressive model of time-varying densities and the necessary theoretical background for understanding the model. Section 3 explains how we may estimate the model and develops the relevant statistical theory. In particular, our estimators of the autoregressive operator and of the error variance operator are shown to be consistent under appropriate regularity conditions. Section 4 presents density forecast along with its asymptotic theory, and provides a variance decomposition method that is useful in studying the moment dynamics of state distributions. In Section 5, we apply our model and methodology to the state distributions representing intra-month distributions of high frequency GBP/USD exchange rate returns and cross-sectional distributions of the NYSE stocks monthly returns. Section 6 summarizes our simulation results on the finite sample performance of the density predictor. Section 7 concludes the paper. All mathematical proofs are collected in Appendix.
In this section, we present our model and some preliminaries that are necessary for the development of our theory and methodology in the paper.
In this paper, we consider a random sequence $(f_t)$ of probability densities of state distributions, which represent some distributional aspect of an economy, for time $t=1,2,\ldots$. The density functions $(f_t)$ are time varying, and are regarded as random elements taking values in the Hilbert space $L^2(C)$ of square integrable functions on some compact subset $C$ of $\mathbb{R}$ or $\mathbb R$ itself. We may introduce the underlying probability space $(\Omega, \mathcal{F}, \mathbb{P})$, and more formally define $f_t:\Omega\to L^2(C)$ for $t=1,2,\ldots$ with the required measurability condition.
We assume that $(f_t)$ has a well defined common mean $\mathbb{E} f$, whose precise meaning will be introduced later in Section (ref). Moreover, we let the demeaned state densities $(w_t)$, $w_t= f_t - \mathbb{E} f$, be generated according to
for $t = 1, 2, \ldots$, where $A$ is a bounded linear operator, called the autoregressive operator, and $(\varepsilon_t)$ is a functional white noise process.\footnote{Of course, we may consider the model $f_t = \mu + Af_{t-1} + \varepsilon_t$ given directly in terms of the state densities $(f_t)$ themselves instead of the demeaned state densities $(w_t)$, where $\mu$ is an additional unknown parameter function corresponding to a constant term in the usual regression. In this paper, we use a more compact presentation in (ref) to focus on the estimation and interpretation of the operator $A$.} The meaning of a functional white noise process will be introduced later in Section (ref). Note that $(w_t)$ takes values in a subset $H$ of $L^2(C)$, which is given by\footnote{Here and elsewhere in this paper, we use the notation $1$ to denote the constant function taking value $1$, instead of the identity function.} \[ H = \left\{v\in L^2(C)\big|\langle 1, v\rangle = 0\right\}, \] where the inner product in the Hilbert space $H$ is defined by $\langle z, w\rangle = \int_C z(x)w(x)dx$, and $1$ denotes the constant function taking value $1$ everywhere on the support. Throughout the paper, we assume that $A$ is an operator on $H$, so that $(\varepsilon_t)$ also take values in $H$.
The model introduced in (ref) may be regarded as a special example of functional autoregression (of order 1), which is often denoted by FAR(1). From this view, our time series $(w_t)$ of demeaned state densities is just a special functional autoregressive process (of order 1). We may define more general FAR($p$) models and processes in a similar way. FAR has the same structure as, and may thus be regarded as, a generalization of the vector autoregression (VAR) that has been widely used in time series econometric modeling. In our FAR model, we simply have functional variables, which are generally infinite dimensional. This contrasts to VAR, where we have a finite number of variables included in the model. FAR and VAR are the same in many aspects such as motivation and intended interpretation. However, technically it is more difficult and involved to deal with FAR than with VAR. For instance, the infinite dimensionality of FAR introduces the ill-posed inverse problem, which makes it harder to do inference. Interested readers are referred to bosq-00 and the references therein for more detailed exposition on FAR.
One advantage of our FAR model for state densities introduced in (ref) is that it may be used to investigate the intertemporal dynamics of the moments of state distributions represented by their densities. To explain our approach, let $v\in H$, and consider a coordinate-process version of our model given as
where $A^*$ is the adjoint of $A$ and $\big(\varepsilon_t(v)\big)$ is a scalar white noise process. The joint operator of a linear operator $A$ on a Hilbert space $H$ is the linear operator $A^*$ on $H$ such that for any $z, w\in H$, $\langle z, Aw\rangle = \langle A^*z, w\rangle$. For any $v\in H$ given, we may interpret $A^\ast v$ as the response function of $\langle v,w_t\rangle$ to an impulse to $w_{t-1}$ given by a Dirac delta function. Note that $\langle v,w_t\rangle$ would increase by $\langle A^\ast v,\delta_x\rangle = (A^\ast v)(x)$, in response to an impulse to $w_{t-1}$ given by the Dirac delta function $\delta_x$ with a spike at $x$.
We may use (ref) to analyze the intertemporal dynamics of various moments of state distributions.\footnote{Here by convention we extend the inner product in $H$ to define $\langle v,w_t\rangle$ for $v\in L^2(C)$. This convention will be made throughout the paper.} For $\iota_p$ defined as $\iota_p(x) = x^p$ for $p=1,2,\ldots$, we have
which is the $p$-th moment of the demeaned state distribution at time $t$. Furthermore, if we write $A^*\iota_p = \sum_{q=1}^\infty c_{p,q}\iota_q$ for a given $p$ with some real sequence $(c_{p,q})$ for $q=1,2,\ldots$, then it follows that
which is an infinite linear combination of all integral moments of the lagged demeaned state distribution represented by $(w_{t-1})$. Consequently, it follows from (ref) that \[ \int_C x^p w_t(x)dx = \sum_{q=1}^\infty c_{p,q} \int_Cx^qw_{t-1}(x)dx + \varepsilon_t(\iota_p), \] which can be used to analyze the moment dynamics of state distributions.
Two special cases appear to be worth mentioning. In what follows, we let $(\mathcal F_t)$ be a filtration, to which $(w_t)$ is adapted. If we assume $A^\ast\iota_2 = \alpha\iota_2$ for some constant $0<\alpha<1$, it follows that $\langle\iota_2, w_t\rangle = \alpha\langle\iota_2, w_{t-1}\rangle + \varepsilon_t(\iota_2)$, from which we may deduce that \[ \mathbb E\left(\langle\iota_2, w_t\rangle\big|\mathcal F_{t-1}\right) = \alpha\langle\iota_2, w_{t-1}\rangle. \] Therefore, (ref) with $v=\iota_2$ reduces essentially to an ARCH model. If in addition to $A^\ast\iota_2 = \alpha\iota_2$, we let $A^\ast\iota_1 = \alpha_1\iota_1 + \alpha_2\iota_2$, then we have \[ \langle\iota_1,w_t\rangle = \beta_1\langle\iota_1,w_{t-1}\rangle + \beta_2\mathbb E\left(\langle\iota_2,w_t\rangle\big|\mathcal F_{t-1}\right) + \varepsilon_t(\iota_1) \] where $\beta_1 = \alpha_1$ and $\beta_2 = \alpha_2/\alpha$. In this case, (ref) with $v=\iota_1$ becomes an ARCH-M model studied in engle-lilien-robins-87. Clearly, our model provides a much more general framework within which we analyze various moment dynamics of general state distributions.
In addition, we may have \[v(x) = 1_B(x),\] where $1_B$ is an indicator function for a subset $B$ of $C$. In this case, $\langle v, w_t\rangle = \int_B w_t(x)dx$ in (ref) represents the probability of the underlying state being at $B$, in terms of the deviation from its expected value. More specifically, for instance, the probability of the underlying state being at a left tail can be analyzed if we set $B=(-\infty,-\tau)\cap C$ with $\tau>0$. We can also choose $B=\big((-\infty,-\tau)\cup (\tau,\infty)\big)\cap C$, so that the probability of the underlying state being at an extreme value can be studied.
Throughout this paper we make the following assumptions:
We use $\left\lVert\cdot\right\rVert$ to denote the operator norm for bounded operators on $H$. The operator norm of $A$ on $H$ is defined as $\left\lVertA\right\rVert = \sup_{v\in H} \left\lVertAv\right\rVert/\left\lVertv\right\rVert$. Recall that a linear operator is compact if it maps the open unit ball to a set that has a compact closure. Since a compact operator on a Hilbert space can be represented as the limit (in operator norm) of a sequence of finite dimensional operators, many features of the finite dimensional matrix theory naturally generalize to compact linear operators. As a consequence, we may regard $A$ essentially as an infinite dimensional matrix. As shown in bosq-00, $\|A^k\|<1$ for some $k\ge 1$ if and only if $\|A^k\|\leq ab^k$ with some $a>0$ and $0<b<1$ for all $k\ge 0$. This condition implies that FAR($1$) introduced in (ref) has a unique stationary solution for $(w_t)$. The i.i.d. assumption for $(\varepsilon_t)$ introduced above is more restrictive than is necessary and can be relaxed to less stringent conditions. Many of our subsequent results hold for sequences that are only serially uncorrelated.
For an $H$-valued random variable $w$, we define its mean by the element $\mathbb{E} w\in H$ such that we have, for all $v\in H$,
If $\mathbb{E}\|w\|<\infty$, then $\mathbb{E} w$ exists and is unique. It can be shown that many of the properties of the usual expectation hold in the functional random variable case. For example, the expectation operator is linear. In particular, if $A$ is a bounded linear operator in $H$, we have that $\mathbb{E}(Aw) = A\mathbb{E} w$. In addition, we have that $\left\lVert\mathbb{E} w\right\rVert \leq \mathbb{E} \left\lVertw\right\rVert$. We define the covariance operator of two $H$-valued random variables $z$ and $w$ by the linear operator $\mathbb{E}(z\otimes w)$ on $H\times H$, for which we have
for all $v\in H$. Note that we have utilized the tensor product notation defined by $(u\otimes v)(\cdot) = \langle v,\cdot\rangle u$ for any $u, v\in H$. If $\mathbb{E}\left\lVertz\right\rVert^2<\infty$ and $\mathbb{E}\left\lVertw\right\rVert^2<\infty$, then $\mathbb{E}(z\otimes w)$ exists and is unique. Naturally, we may call $\mathbb{E}(w\otimes w)$ the variance operator of $w$.
An $H$-valued random variable $w$ is called Gaussian if the random variable $\langle v, w\rangle$ is Gaussian for all $v\in H$. We denote by $\mathbb{N}(\mu, \Sigma)$ a Gaussian random variable with mean $\mu$ and variance operator $\Sigma$.
An operator $A$ on a separable Hilbert space $H$ is compact if and only if it can be written as
for some orthonormal systems $(u_k)$ and $(v_k)$ of $H$ and a sequence $(\kappa_k)$ of nonnegative numbers tending to 0. The compact linear operator is called nuclear if $\sum_{k=1}^\infty \left\lvert\kappa_k\right\rvert < \infty$, and Hilbert-Schmidt if $\sum_{k=1}^\infty \kappa_k^2 < \infty$. It is well known that $\mathbb{E}(z\otimes w)$ is nuclear, and therefore compact. Furthermore, the variance operator $\mathbb{E}(w\otimes w)$ is self-adjoint and nonnegative (strictly positive if non-degenerate), and admits the spectral representation
for some nonnegative (strictly positive if non-degenerate) sequence $(\lambda_k)$ and some orthonormal basis $(v_k)$ of $H$.
The spectral representation of the autoregressive operator as in (ref) gives rise to interesting interpretations of the autoregressive operator. Note that, under (ref), our model in (ref) can be written as
If we order $(\kappa_k)$ and correspondingly $(v_k)$ and $(u_k)$ such that $\kappa_1>\kappa_2>\ldots$, then $v_1$ may be viewed as the direction in which $w_{t-1}$ mainly affects $w_t$. The corresponding feature of the distribution generated from $v_1$ will be called the leading progressive feature. On the other hand, $u_1$ represents the direction in which $w_t$ is mainly affected by $w_{t-1}$. The corresponding feature of the distribution generated from $u_1$ will be called the leading regressive feature. For example, if $v_1 = \iota_2$ and $u_1 = \iota_1$ so that the leading progressive feature is the second moment and the leading regressive feature is the first moment, then in this process the second moment of $w_{t-1}$ mainly affects the distribution of $w_t$ and the first moment of $w_t$ is mainly affected by the distribution of $w_{t-1}$.
A sequence $(\varepsilon_t)$ of $H$-valued random variables is called a white noise if $\mathbb{E}(\varepsilon_t\otimes \varepsilon_t)$ is the same for all $t$ and $\mathbb{E}(\varepsilon_t \otimes \varepsilon_s)=0$ for all $t\neq s$. It is easy to see that if $(\varepsilon_t)$ is a white noise, then for any $v\in H$, the scalar process $(\varepsilon_t(v))$, $\varepsilon_t(v) = \langle v, \varepsilon_t\rangle$, is a white noise. Under Assumption (ref), the model (ref) has a unique stationary solution for $(w_t)$, which is given by
where $(\varepsilon_t)$ is a white noise process. See, for instance, bosq-00 for more details. It is therefore convenient to define the autocovariance function of $(w_t)$ by
for $k\in \mathbb Z$, as in the analysis of finite dimensional time series. Of particular interest are $P = \Gamma(1)$ and $Q = \Gamma(0)$. It is easy to deduce that $\Gamma(k) = A^kQ$ for $k\geq 1$ and $\Gamma(k) = \Gamma(-k)^*$ for all $k$. In particular, we have
which will be used to estimate $A$.
It is well known that compact linear operators on infinite dimensional Hilbert spaces are not invertible. Therefore one may not directly use the relationship in (ref) to define the autoregressive operator $A$ as $A = PQ^{-1}$ since $Q^{-1}$ is not well defined on $H$. In fact, if the kernel of $Q$ is $\{0\}$, then $Q^{-1}$ is well defined only on $\mathcal{R}(Q) = \{u\in H|\sum_{k=1}^\infty \langle u, v_k\rangle^2/\lambda_k^2 <\infty\}$, which is a proper subspace of $H$. Consequently, we have $A=PQ^{-1}$ well defined only on the restricted domain $\mathcal{D}=\mathcal{R}(Q)$. This problem is often referred to as the ill-posed inverse problem.
The standard method to circumvent this problem in the functional data analysis literature is to restrict the definition of $A$ in a finite dimensional subspace of $H$. To explain this method, let $\lambda_1>\lambda_2>\cdots>0$ and define $V_K$ to be the subspace of $H$ spanned by the $K$ eigenvectors $v_1, \ldots, v_K$ associated with the eigenvalues $\lambda_1, \ldots, \lambda_K$ of $Q$. Let $Q_K = \Pi_KQ\Pi_K$ where $\Pi_K$ denotes the projection onto $V_K$ and define
i.e., the inverse of $Q$ on $V_K$.\footnote{See, e.g., benatia-carrasco-florens-17 for the use of other methods of regularization.}
Now let
which is the autoregressive operator $A$ restricted to the subspace $V_K$ of $H$. Note that $V_K$ is generated by the first $K$ principal components of $(w_t)$. Since $(\lambda_k)$ decreases to zero, we may expect that $A_K$ approximates $A$ well if the dimension of $V_K$ increases. The estimator of $A$, which will be introduced in the next section, is indeed the sample analogue estimator of $A_K$ in (ref), and we let $K$ increase as $T$ increases.
In most applications, $(f_t)$ are not directly observed and should therefore be estimated before we look at the FAR model specified by (ref). We suppose that $N$ observations from the probability density $f_t$ are available so that we may estimate $f_t$ consistently for each $t = 1, 2, \ldots$. In what follows, we denote by $\hat f_t$ the consistent estimator of $f_t$ based on $N$ observations and define $\Delta_t = \hat f_t - f_t$ for $t = 1, \ldots, T$. As before, we let $(\lambda_k)$ be the ordered eigenvalues of $Q$.
Assumption (ref) is very mild and appears to be widely satisfied in practical applications. The condition in part (ref) holds if and only if the kernel of $Q$ contains only the origin, which is met as long as the distribution of $(f_t)$ is non-degenerate. Boundedness of $(f_t)$ in part (ref) is expected to hold in many practical applications, since $(f_t)$ is a family of density functions. This is not strictly necessary, and is only required for an optimal rate in the strong laws of large numbers we establish in Section 3.1 and 3.2. For unbounded $(f_t)$, we have a reduced rate of convergence for the strong laws of large numbers. However, the $L^2$-convergence results in Section 3.1 and 3.2 hold for unbounded $(f_t)$ without any modification. The condition in part (c) holds in many applications. For instance, if $(f_t)$ are estimated for each $t$ by the kernel density estimator from $N$ observations, we may expect for all $t$ that $\mathbb{E}\left\lVert\Delta_t\right\rVert^2 = O(N^{-r})$ for some $0<r<2/5$ under very general regularity conditions which allow in particular for dependency among $N$ observations. See, e.g., bosq-98. One may expect even higher rates of convergence if higher-order kernels are used. See, e.g., wand-jones-95 .
First we estimate $\mathbb{E} f$ by the sample average of $(f_t)$, i.e., by
The following result holds.
Lemma (ref) establishes the $L^2$ and a.s. consistency for the sample mean $\bar f$. It shows in particular that using estimated densities $(\hat f_t)$ in place of $(f_t)$ does not affect the convergence rates established in Theorem 3.7 and Corollary 3.2 of bosq-00 as long as the number of observations used to estimate $(f_t)$ is sufficiently large. The conditions in Assumptions (ref) and (ref) are sufficient to yield optimal rates in the laws of large numbers here.
Similarly, we may use the sample analogue to consistently estimate the autocovariance operators $Q$ and $P$. Let $\hat w_t = \hat f_t - \bar f$ for $t = 1, 2, \ldots$, and define
and
Then we have
Theorem (ref) shows that the sample analogue estimators $\hat Q$ and $\hat P$ based on the estimated density functions $(\hat f_t)$ are consistent, and have the same rates of convergence as those based on the true density functions $(f_t)$, which have been established in Theorem 4.1 and Corollary 4.1 of bosq-00.
The implementation of our methodology requires the estimation of eigenvalues and eigenvectors of $Q$, and $Q$ can be consistently estimated by $\hat Q$ by Theorem (ref). Naturally, the eigenvalues and eigenvectors of $Q$ are estimated by those of $\hat Q$, which we denote by $(\hat \lambda_k, \hat v_k)$. We assume that the estimated eigenvalues are distinct and order them so that $\hat \lambda_1> \hat \lambda_2 > \cdots$\footnote{We could potentially allow for multiplicity. However, in that case the eigenvectors could not be uniquely identified even after normalization. As a consequence we need to show consistency of eigenspaces, which complicates the notations and proofs substantially. In this paper, we assume for convenience that the eigenvalues are distinct.}. The eigenvalues and eigenvectors of $\hat Q$ are the pairs $(\hat \lambda, \hat v)$ such that $\hat Q\hat v = \hat \lambda \hat v$.
We may solve the eigen-problem either by discretizing $(\hat w_t)$ or by expanding $(\hat w_t)$ into linear combinations of an orthonormal basis of $H$ and obtaining the matrix representation of $\hat Q$ with respect to that basis. Either way, we transform the original problem into the problem of computing eigenvalues and eigenvectors of a matrix. Readers are referred to Section 8.4 in ramsey-silverman-05 for more details.
For definiteness, let us assume that the eigenspace corresponding to $(\lambda_k)$ is one dimensional for each $k$ from now on. This implies that Assumption (ref)(ref) should be replaced by $\lambda_1>\lambda_2>\ldots>0$. We define $v_k'=\operatorname*{sgn}\langle \hat v_k, v_k\rangle v_k$. Note that since both $(v_k)$ and $(-v_k)$ are eigenvectors corresponding to $\lambda_k$, the introduction of $(v_k')$ is essential for definitiveness of eigenvectors. The following lemma from bosq-00 associates $\left\lvert\hat \lambda_k - \lambda_k\right\rvert$ and $\left\lVert\hat v_k - v_k'\right\rVert$ with $\left\lVert\hat Q - Q\right\rVert$.
Let $\tau(K) = \sup_{1\leq k \leq K}(\lambda_k - \lambda_{k+1})^{-1}$ and set $K$ to be a function of $T$, i.e., $K_T$ such that $K_T\to \infty$ as $T\to \infty$. It follows directly from Theorem (ref) and Lemma (ref) that
Our results in Theorem (ref) shows that the convergence rates of the estimated eigenvalues, in the presence of estimation errors for the densities, are the same as in Theorem 4.4 and Corollary 4.3 in bosq-00.
Using the estimators introduced in Sections (ref) and (ref), we may consistently estimate the autoregressive operator $A$ by an estimator of $A_K$ given earlier in (ref). To be explicit, we define
and subsequently
Recall that $K$ is a function of the sample size $T$ such that $K\to \infty$ as $T\to \infty$. Theorem (ref) establishes the consistency of $\hat A_K$. Consistency as in Theorem 8.8 in bosq-00 continues to hold in the presence of estimation errors for the densities.
Now we define
where $(\hat \varepsilon_t)$ obtained by
are the fitted residuals from the FAR in (ref) with $(w_t)$ replaced by $(\hat w_t)$. As is expected from Theorem (ref), we have
Corollary (ref) establishes the strong consistency of $\hat \Sigma$.
This section provides some tools and asymptotics to analyze the distributional autoregression. In particular, we discuss two topics including the forecast of state density and the moment dynamics of state distributions.
Our model can be used to obtain the forecasts of future state densities. For one-step forecast, we use
where $\hat A_K$ is the estimator of the autoregressive operator $A$ introduced in (ref). Multiple step forecasts can be obtained by recursively applying $\hat A_K$ to the forecasted values in the previous steps. Below we develop the asymptotic theory for our one-step forecast.
The development of our asymptotics require some extra technical conditions in addition to what we have in Assumptions (ref) and (ref).
Condition (ref) is weak and satisfied by many sequences including $\lambda_k = c/k^a$, $\lambda_k = c/(k^a\ln^b k)$ and $\lambda_k = ce^{-ak}$, where $a, b$ and $c$ are some positive constants. Let $\delta_k = \lambda_k - \lambda_{k+1}$. Then $\delta_k$ is a decreasing sequence by convexity of $\lambda_k$. The condition (ref) is met whenever the tail probability of $\langle v_k, w_t\rangle$ decreases fast enough. For instance, if $w_t$ is Gaussian, the condition holds with $M=3$. Note that $\langle v_k, w_t\rangle$ has mean zero and variance $\lambda_k$ for $k = 1, 2, \ldots$. It may not hold if the distribution of $\langle v_k, w_t\rangle$ for some $k$ has a thicker tail.
In the development of our subsequent theory, we denote by $\Pi_K$ the projection onto the subspace spanned by the first $K$ principal eigenvectors $v_1, \ldots, v_K$ of $Q$, and denote by $\hat \Pi_k$ the projection onto the subspace spanned by the first $K$ principal eigenvectors $\hat v_1, \ldots, \hat v_K$ of $\hat Q$. The asymptotic theory that we develop for the one-step forecast is based on the foundational results due to mas-07. To avoid technical difficulties, we follow mas-07 to compute the estimator $\hat A_K$ using samples only up to time $T-1$.
The condition (ref) requires that $K$ should not increase too fast relative to $T$. Loosely put, condition (ref) requires that $A$ should be at least as smooth as $Q^{1/2}$. If indeed they have common eigenvectors, then this condition holds if and only if the ordered sequence of eigenvalues of $A$ decays at least as fast as the ordered sequence of eigenvalues of $Q^{1/2}$.
The result in Lemma (ref) is not directly applicable to obtain the asymptotic confidence interval for the forecast of $(w_t)$. In particular, it has a random bias term $(A\hat \Pi_K -A) w_T$. This can be remedied, at the expense of a set of additional conditions as in the following theorem.
The newly introduced conditions in Theorem (ref) can be restrictive, and depend upon the decaying rate of the eigenvalues of the variance operator $Q$, which is unobserved, and asymptotic independence of the process $(w_t\otimes \Delta_t)$. However, they are not prohibitively stringent. For instance, the conditions on the eigenvalues are satisfied for $\lambda_k = ce^{-k}$ where $c$ is some positive constant, and it seems that in many applications the eigenvalues indeed exhibit an exponential decay. The conditions on the correlations are satisfied if $(w_t\otimes \Delta_t)$ is a mixingale or $\rho$-mixing, with geometrically decaying mixing rates.
The asymptotic result in Theorem (ref) can be used to construct the asymptotic confidence interval for the one-step forecast of $(w_t)$. Indeed, we may easily deduce from the result in Theorem (ref) that
It follows that for any $v\in H$ the interval forecast for $\langle v, w_{T+1}\rangle$ with confidence level $\alpha$ is given by
with $z_{\alpha/2}=\Phi^{-1}(1-\alpha/2)$, where $\Phi$ is the cumulative distribution function of the standard normal. Here $\Sigma$ could be replaced by its consistent estimate in implementation.
In this section, we demonstrate how we may analyze the moment dynamics of state distributions. The analysis proves useful in learning the dependence structure in the moments of state distributions across time. For the analysis, we define a normalized moment basis $(u_k)$ in $H$ such that $u_k$ is the $k$-th order polynomial for $k=1,2,\ldots$, and that
Such a basis could be obtained by the Gram-Schmidt orthogonalization procedure. For $C = [a,b]$, we may easily see that $u_1$ is given by $u_1(x) = C_1\big[x-(a+b)/2\big]$, where $C_1$ is a constant determined by $\langle u_1, Qu_1\rangle = 1$. Moreover, once $u_1, u_2, \ldots, u_p$ are obtained, we may readily find $u_{p+1}$ such that $\langle 1,u_{p+1}\rangle = 0$, $\langle u_q,Qu_{p+1} \rangle = 0$ for $q=1,\ldots,p$ and $\langle u_{p+1}, Q u_{p+1}\rangle = 1$. Note that $(u_k)$ is orthonormal with respect to the inner product $\langle \cdot, Q\cdot\rangle$ and generates $H$. That is, it is an orthonormal basis of $H$. Since $u_k$ is a polynomial of order $k$, we shall view it as the function that generates the normalized $k$-th order moment of state distributions.
Now we consider (ref), which defines the dynamics for the $v$-moment of the underlying state distributions. For any $v\in H$ given, write \[ A^\ast v = \sum_{k=1}^\infty \langle A^\ast v,Qu_k\rangle u_k = \sum_{k=1}^\infty \langle v,AQu_k\rangle u_k, \] so that we may deduce
which are quite useful to analyze the moment dynamics of state distributions. Note that we have \[ \mathbb E\langle u_p,w_t\rangle\langle u_q,w_t\rangle = \langle u_p,Qu_q\rangle = \left\{
\right. \] by construction. Consequently, for any $v\in H$ given, we have
due to (ref) and (ref).
In our analysis of the moment dynamics of state distributions, for any given $v\in H$ we define
to be the $R^2$ of the $v$-moment of state distribution, and interpret the ratio
as the proportion of the variance of the $v$-moment of the current state distribution contributed by the $k$-th moment of the past state distribution. Clearly, $R_v^2$ and $\pi_v(k)$ defined respectively in (ref) and (ref) can be consistently estimated by \[ \hat R_v^2 = 1-\frac{\langle v,\hat\Sigma v\rangle}{\langle v,\hat Qv\rangle} \quad\mbox{and}\quad \hat\pi_v(k) = \frac{\langle v,\hat A_K\hat Qu_k\rangle^2}{\langle v,\hat Qv\rangle} \] for any $v\in H$, where $\hat\Sigma, \hat Q$ and $\hat A_K$ are the estimates of $\Sigma$, $Q$ and $A$, respectively, which have been introduced in Section 3.
In this section we use our methodology to analyze two financial markets. In the first, we study the intra-month distributions of the GBP/USD exchange rate 15-minute log returns. In the second, we study the cross-sectional distributions of monthly returns of stocks listed on the NYSE. We explore various aspects of the distributional dynamics in the foreign exchange and stock markets that would not be revealed if we only look at their aggregate time series.
We also present the forecast performances of our model. For each forecast period, we obtain the predicted density $\hat f$ and calculate how much it deviates from the actually observed density $f$, using the six measures given in Table (ref). $\hat F$ and $F$ in the table are the distribution functions corresponding to $\hat f$ and $f$, respectively. The first two measures, namely the $L^2$- and the $L^1$-deviations, are based on the density functions. The middle two measures, namely the Kolmogorov-Smirnov and the Cram\'{e}r-von Mises deviations, are based on the distribution functions. The last two measures, namely the (absolute) deviation in mean and in variance, are based on the first two moments of the densities. Notice that all these quantities are always nonnegative, and are zero when $\hat f = f$.
For a high frequency intra-month financial return series, we may consider the individual intra-month returns as draws from a common distribution, which varies by month, under the piecewise stationarity assumption. These draws may be dependent, hence accommodating the case in which the intra-month returns are not independent but are strictly stationary. This view not only provides a justification for applying our model to time-varying return distributions, but also reflects the fact that the intra-month returns exhibit statistical characteristics that do not carry over to longer horizons.
In this application, we look at the intra-month distributions of the GBP/USD exchange rate 15-minute log returns. Since the foreign exchange market operates globally around-the-clock on weekdays, it is more natural to treat every four weeks as a month than to follow the calendar months. We thus set a month in this application to be a four week period. We use data from January 4th, 1999 to April 3rd, 2015, and split the data into a total of 212 months using the four-week convention. The number of observations in each period ranges from 1550 to 1904, with a mean of 1880. We obtain the kernel density estimate for each period using the Epanechnikov kernel with the optimal bandwidth $h_t=2.3449\hat \sigma_t n_t^{-1/5}$, where $\hat \sigma_t$ and $n_t$ are respectively the sample standard deviation and the sample size in period $t$. The support of the densities is set to be $[-0.0043, 0.0043]$, which includes more than $99.9\%$ of the actual observations. We plot the time series of both the densities and the demeaned densities in Figure (ref), from which we could easily see the time-varying nature of the corresponding distributions.
To implement our method, we represent the densities with the Daubechies wavelets using 1037 basis functions. Each density is then represented as a 1037-dimensional vector whose coordinates are the wavelet coefficients of the density function. Thus we transform the functional time series into a time series of high dimensional vectors. We obtain the matrix representations of the sample variance operator $\hat Q$ and the sample first-order autocovariance operator $\hat P$ with respect to that basis, which are respectively the sample covariance and the sample first-order autocovariance matrices of the vector time series. We then estimate the eigenvalues and eigenvectors of the variance operator $Q$ respectively by the eigenvalues and eigenvectors of the sample covariance matrix of the vector time series. When we approximate the inverse of the variance operator, we set $K=4$ to obtain the best rolling out-of-sample forecast performance (see below). This is the data-driven method we propose for choosing the value of the nuisance parameter $K$. As another justification for our choice of $K$ from the perspective of functional principal component analysis, we note that the first 4 principal components together explain 99.7% of the total variation in the density process. See Figure (ref) for the scree plot of the ten largest eigenvalues of the variance operator.
As an illustration of how we may use the tools developed in this paper to perform interesting analyses, we first obtain the leading progressive and regressive features for this density process. Each feature is represented by its values at 1024 points evenly spaced on the support of the distributions. See Figure (ref). The leading progressive feature is a concentrating feature: If it is loaded with a positive scale, the distribution will be more concentrated around its center. It also shows that the relative frequency of observations with small returns is most important for determining distributions in the future. The leading regressive feature shows it is also the relative frequency of observations with small returns that is most affected by distributions in the past. Both features indicate that the observations with extreme returns are not important in the distributional dynamics.
We also illustrate how various aspects in the current distribution respond to shocks to the past distribution. The left two panels in Figure (ref) present the current first two moments' responses with respect to Dirac-$\delta$ impulses to the previous distribution at different levels of return. The shaded areas give the corresponding 95% residual bootstrap confidence bands based on 2000 repetitions. The top left panel shows that the mean of the current distribution is most affected by shocks to the relative frequencies of the observations with small-to-moderate returns in the previous period, and shocks to the relative frequencies of the observations with large or extreme returns in the previous period are not important. Moreover, we see the momentum effect, i.e., a larger number of observations with small-to-moderate negative (positive) returns in the previous period is likely to result in a lower (higher) first moment in the current period. The bottom left panel shows that the second moment of the current distribution, on the other hand, is affected by shocks to the relative frequencies of returns at all levels in the previous period. One may expect a smaller second moment in the current period if there were more observations with small-to-moderate returns in the previous period, and a larger second moment in the current period if there were more observations with large returns in the previous period. The effects of small-to-moderate returns in the previous period on the second moment in the current period are not monotonic. The stabilizing effects of some moderate returns are at least as strong as those of the small returns.
The right two panels give the proportions of the variances of the current first and second moments that are explained by the variances of the first ten moments of the previous distribution, with their 95% residual bootstrap confidence intervals. The corresponding $R_v^2$ are respectively $0.0380$ and $0.8091$, with 95% bootstrap confidence intervals $[0.0085, 0.1186]$ and $[0.6214, 0.8727]$. Though many integral moments, especially those of odd orders, affect the current mean, their overall effect is almost negligible. In contrast, the past second moment, whose variance explains about $80\%$ of the total variance, is very informative about the current second moment, while the other past moments provide almost no information. Our findings here are entirely consistent with the widespread and repeated observations by many empirical researchers that in financial markets the mean is usually difficult to predict, while the volatility process is much more persistent and therefore much more predictable.
Figure (ref) gives the impulse responses and the variance decompositions of the tail probabilities. The left tail includes the returns that are smaller than the 5th percentile of the mean density, and the right tail includes the returns that are greater than the 95th percentile of the mean density. In this application, the responses of the two tail probabilities to shocks in the past density are quite symmetric. Shocks to, and subsequent changes in, the relative frequencies of returns at all levels in the previous period matter. A larger number of observations with small-to-moderate returns in the last period tends to presage thinner tails of the current distribution, while a larger number of observations with large returns in the previous period tends to presage thicker tails of the current distribution. The tail responses of small-to-moderate returns are non-monotonic, and are essentially the same as their responses to the second moment. In the variance decomposition analysis, about $80\%$ of the variances of the current tail probabilities are explained by the variances of the past second moment, which shows that the past second moment is very informative about the current tail probabilities.
Finally, to evaluate the prediction performance of our model in this application, we make rolling out-of-sample forecasts of the state densities for the last 50 periods. For each forecast period, we use the data prior to that period to estimate the autoregressive coefficient. We dynamically choose $K$ by cross-validation, using the five periods before the forecast period as the validation periods. We then make one-period-ahead forecasts of the demeaned density as in (ref). We add the mean density to the forecasted demeaned density to obtain the density forecast. We normalize the forecasted density by setting its negative values to zero and then rescaling it so that it integrates to one. In the end, we calculate the forecast error. For comparison, we consider two other predictors, which are respectively the time average of all the estimated densities and the estimated density of the last period. We label our functional autoregressive predictor as the FAR predictor, and the other two respectively as the AVE predictor and the LAST predictor. The AVE and the LAST predictors are used as benchmarks: If the true process of the random densities is indeed independent and identically distributed over time, one would expect the AVE predictor to perform well. If the true process is a nonstationary martingale process, one would expect the LAST predictor to perform the best. Comparisons with the two benchmarks reveal information on how far away the actual process is from the i.i.d. scenario and from the nonstationarity scenario.
Table (ref) reports the means and medians (in parentheses) of the 50 periods of forecast errors for each of the three predictors. In terms of the first four measures, the FAR predictor performs roughly the same as the LAST predictor, and they outperform the AVE estimator significantly. This suggests that the true process does not resemble an i.i.d. process, and that there is likely to be strong persistency in the density process. In terms of predicting the mean, the FAR predictor performs about the same as the AVE predictor, and much better than the LAST predictor. In terms of predicting the variance, the FAR predictor performs about the same as the LAST predictor, and much better than the AVE estimator. These results corroborate the findings that in financial markets the mean is usually difficult to predict, while the volatility is quite persistent, and therefore much more predictable.
In this application, we look at the cross-sectional distributions of monthly returns of stocks listed on the New York Stock Exchange. In each month, the individual stock returns are viewed as draws from a common distribution. This distribution changes from month to month and is assumed to follow our FAR model (ref). We use data from January of 1980 to December of 2014, a total of 420 months. In each month, we only take into account the stocks that are traded. The number of observations each month ranges from 1926 to 3076, with a mean of 2464. We use the same method as in the previous application to estimate the densities. The support of the density in this application is set to be $[-0.6071, 1.0548]$, which includes more than $99.9\%$ of the actual observations. Figure (ref) plots the times series of the densities and the demeaned densities of the NYSE stocks monthly returns.
We use the same procedures as in the previous application for model estimation. We set $K=3$ based on the out-of-sample forecast performance. Figure (ref) gives the scree plot of the eigenvalues of the sample variance operator in this application. The first 3 principal components explain $97.0\%$ of the total variation in the density process. Figure (ref) demonstrates the leading progressive and regressive features for the density process of the NYSE stocks monthly returns. As in the GBP/USD exchange rate application, it is the relative frequency of small returns that is most critical in the distribution dynamics. The small returns play a central role in affecting the future distribution, and they are most affected by the distribution in the previous period.
In Figure (ref), we provide the impulse responses and the variance decompositions of the first two moments in this application. The corresponding $R_v^2$ are $0.0368$ and $0.3376$ respectively with 95% residual bootstrap confidence intervals $[0.0108, 0.0815]$ and $[0.2199, 0.4709]$, which implies that the first moment is much less predictable than the second moment through the integral moments of the past distribution. Unlike the exchange rate returns application, in which the current mean responds only to changes in the relative frequencies of small returns in the previous period, the current mean in this application responds more significantly to changes in the relative frequencies of moderate-to-large returns in the previous period. The typical momentum effect is generated by stocks with moderately large returns, both positive and negative. For moderately small negative returns, we even find some evidence of a reverse momentum effect, a pattern which is also found in international stock returns by fama-french-12. As expected, the second moment of current return distribution decreases as we have more small return stocks in the past distribution. Interestingly, the current second moment responds positively to changes in the relative frequencies of moderate-to-large negative returns, but not to changes in the relative frequencies of positive returns, in the previous period. This is presumably due to the leverage effect, and our result here suggests that the leverage effect is asymmetric, being effective only for negative return shocks. The past integral moments are not informative about the current mean. They are relatively more informative about the current second moment. However, the relevance is smaller than in the previous application, and the past fourth and first moments appear to play some non-negligible roles.
Figure (ref) presents the impulse responses and the variance decompositions of the tail probabilities in this application. The tail events are defined in the same way as in the last application. In contrast with the previous application, the left tail probabilities and the right tail probabilities in this application respond differently to Dirac-$\delta$ impulses to the last period's density. As is expected, positive shocks to the relative frequencies of small returns decrease both the left and the right tail probabilities in the next period. Positive shocks to the relative frequencies of moderate-to-large negative returns tend to presage a thicker left tail, whereas positive shocks to the relative frequencies of moderate-to-large positive returns tend to presage a thicker right tail. Of the two effects, the former appears to be much more significant than the latter. Moreover, positive shocks to the relative frequencies of moderate-to-large positive returns decrease the left tail probability, whereas positive shocks to the relative frequencies of moderate-to-large negative returns also increase the right tail probability, though less substantially. On the other hand, the current left tail probability is related to both the past first and second moments, while the current right tail probability hinges more on the past second and fourth moments.
In this application, we make rolling out-of-sample forecasts for the last 100 periods since we have longer time series available. Table (ref) summarizes the prediction results. In terms of the first four measures, our FAR predictor performs the best while the LAST estimator performs the worst. In terms of predicting the mean, the FAR predictor performs slightly worse than the AVE predictor, and much better than the LAST predictor. In terms of predicting the variance, the FAR predictor performs worse than the LAST predictor, but much better than the AVE predictor. These results suggest that this density process is likely to be stationary, and our model has good forecast power. However, there may be some particular directions in which the process is close to nonstationary, as volatility appears to be quite persistent.
In this section we study the finite sample properties of our FAR estimator by conducting out-of-sample forecast using simulated data. To make our simulations more realistic and practically relevant, we use the estimation results from the previous two empirical applications and generate data that mimics the estimated empirical models as closely as possible.
To be specific, for each application, we take the estimated coefficient $\hat A_K$ to be the autoregressive coefficient in the FAR(1) data generating process given in (ref) for the process $(w_t)$. In each iteration of simulations, the error process is bootstrapped from the demeaned residuals of the estimated FAR(1) model, given in (ref). Once we obtain the simulated demeaned densities $(w_t)$, we add back the time average of the estimated densities $\bar f$ and normalize them as in the forecast procedures in the two applications to obtain simulated density functions. In each simulation, we start from $\varepsilon_0 = 0$, simulate 1000 periods and take the last $T+1$ periods. The actual values for $T$ will be specified below. We take the time series of the last $T+1$ density functions as the true process of densities. We then simulate $N$ observations from the density in each period using acceptance (or rejection) sampling. The actual values for $N$ will be specified below. Acceptance sampling is a widely used sampling method, especially in Bayesian statistics, when one can not directly sample from the distribution but the corresponding density function is known. The researcher samples from a proposed distribution which is different from the target distribution, and accepts or rejects each draw with a probability that is determined by the ratio of the proposal and target density functions evaluated at the draw. In our sampling procedure, we use the uniform distribution on the support as the proposal distribution. We take the accepted draws as the simulated observations.
We implement our functional autoregression using the simulated observations from the first $T$ periods, leaving the last period for forecast performance evaluation. The estimation procedure as well as the parameter settings are the same as in the empirical applications, except that in the kernel density estimation step, we use the Normal kernel instead of the Epanechnikov kernel, with the corresponding optimal bandwidth. The Normal kernel performs better in our simulations in that it gives smaller kernel density estimation error. Once we obtain the estimated autoregressive coefficient, we do one-step-ahead forecast for the last period and compare our forecast with the true density in that period. The forecast procedure is exactly the same as in the two applications in Section 5. We evaluate the difference of our forecast from the true density using the six criteria introduced in the previous section.
We set $T$ to be 50, 100, 200, or 500, and $N$ to be 100, 200, 500 or 5000. For each pair of $(T, N)$, the mean and median of the forecast errors based on 2000 iterations of simulations are reported in Table (ref) and Table (ref). For convenience, we label the simulations based on the two estimated models as the “FOREX simulations” and the “NYSE simulations” respectively. We make the following observations.
First, for both sets of simulations, the forecast performance is not very sensitive to the number of periods ($T$) we have in the estimation. Once we have at least 100 periods, results become highly stable. It seems that a total of 50 periods is still acceptable, although the forecast errors are slightly bigger.
Second, forecast performances are more sensitive to the number of observations ($N$) in each period. For both sets of simulations and for all values of $T$, when we increase the number of observations in each period, the forecast errors decrease non-negligibly. This suggests that a large part of the forecast error comes from the density estimation procedure. This agrees with the usual observation that kernel density estimation converges at a relatively slow rate.
Third, for both sets of simulations, and for all combinations of $T$ and $N$, the FAR estimator performs best in terms of the first four criteria. This implies that our method has good curve fitting property, in terms of both the density function (the first two criteria) and the distribution function (the middle two criteria). In the FOREX simulations, the FAR predictor performs roughly the same as the AVE predictor, and much better than the LAST estimator in predicting the mean. In predicting the variance, the FAR predictor performs roughly the same as the LAST predictor, and much better than the AVE predictor. In the NYSE simulations, the FAR predictor has the best performance in predicting both the mean and the variance.
To summarize, we find that the FAR predictor forecasts densities accurately, and that the forecast is stable, even with a relatively short time series.
This paper presents a model for time-varying densities that follow an autoregressive stationary process in some function space. We explore the properties and implications of this model, showing how it relates to many familiar time series models such as the ARCH-type models. Our approach is highly nonparametric, in that it does not assume any parametric form for the underlying distributions. We have outlined the procedures to estimate the model assuming that the densities are not observable, and have shown the consistency of the estimator as both the sample size for the density estimation and the time horizon go to infinity. Furthermore, we illustrate how to forecast densities using our model, and establish the asymptotic normality of the predictor. We develop tools such as progressive/regressive analysis, impulse response analysis and variance decomposition to help study the dynamics in various aspects of the distribution process. We discover several interesting features in the distribution process of the intra-month GBP/USD exchange rate 15-minute log returns and in the distribution process of the NYSE stocks monthly returns. Finally, Monte Carlo simulations show that the FAR density predictor outperforms its competitors. In sum, the FAR modeling for time series of state distributions is a very general framework that is highly flexible, readily applicable, and it provides a powerful method to extract information of a distributional process and generate accurate forecasts.