{'fraction_non_alphanumeric': 0.04090909090909091, 'fraction_numerical': 0.010427807486631016, 'mean_word_length': 4.48679706601467, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Testing Of Granger Causality In Functional Time Series

We develop a multivariate functional autoregressive model (MFAR), which captures the cross-correlation among multiple functional time series and thus improves forecast accuracy. We estimate the parameters under the Bayesian dynamic linear models (DLM) framework. In order to capture Granger causality from one FAR series to another we employ Bayes Factor. Motivated by the broad application of functional data in finance, we investigate the causality between the yield curves of two countries. Furthermore, we illustrate a climatology example, examining whether the weather conditions Granger cause pollutant daily levels in a city.


Introduction
A functional time series arises when at each point in time we observe a random function, generally at discrete points on a dense support. Examples of functions include satellite images(Guyet & Nicolas [1]), fMRI data(Stoehr, Aston & Kirch [2]) and yield curves ( Sen & Klüppelberg [3]). Another common practice that leads to functional time series is to slice high-frequency temporal sequence records into curves over logical consecutive time intervals. Examples of such curve time series include the daily price curve of financial transactions ( Horváth, Kokoszka & Rice [4]; Li, Robinson & Shang [5]), electricity price curves as a function of day-of-year ( Chen & Li [6]; Chen, Marron & Zhang [7]), daily pollution levels as a function of the time of the day (Aue, Norinho & Hörmann [8]), and daily patterns of environmental data (Besse, Cardot & Stephenson [9], Ramsay & Dalzell [10]). The functional autoregressive (FAR) model, developed by Bosq [11], plays a central role in modeling and predicting functional time series. The FAR model extends the vector autoregressive (VAR) model to the infinite-dimensional space and considers the serial correlation between the functional observations.
Most of the existing literature relies on functional principal component analysis (FPCA) to treat multivariate functional data. In a multivariate setting, along with serial correlation, one needs to study the correlation among the data series. Ramsay & Silverman [12] introduced a classical multivariate FPCA that concatenates the multiple functional data and performs univariate FPCA on a single long vector. However, this model is restricted to multivariate random functions on the same finite, one-dimensional interval. Chiou, Chen & Yang [13] discussed a normalized multivariate functional principal component analysis (mFPCA n ) which accounts for data measured in different units using normalized covariance operators. Although mFPCA n produces sound estimates on dense grids without calculation errors, it leads to inadmissible results with significant measurement errors in the case of sparse settings.
Apart from FPCA, there is abundant literature that models multivariate functional data using regression analysis. Zhu, Li & Kong [14] developed a multivariate varying coefficients model (MVCM) to estimate covariance function. They also establish the uniform convergence rate of this covariance function. Li et al. [15] extended the previous varying-coefficient single index model to perform the regression on functional response data. Zhu et al. [16] performed linear transformation and adopted the Bayesian approach to sample the model parameters in the transformed space. Zhu, Strawn, & Dunson [17] introduced a notion of graphical models in multivariate functional data by constructing Bayesian inference of undirected, decomposable graphs based on Markov distribution. The literature mentioned above mainly focuses on dense functional data. Li, Xiao & Luo [18] adopted the fast covariance-based FPCA method for multivariate sparse functional data, approximating covariance functions using Tensor-product B-splines. The underlying assumption while using FPCA is that the observable process with measurement errors follows the multivariate functional autoregressive (MFAR) process. We later prove that such an observable process follows a multivariate functional autoregressive moving average (FARMA) process. Thus, existing approaches based on FPCA produce inefficient estimates. This paper proposes a hierarchical model for multivariate functional time series, which is the extension of Kowal, Matteson & Ruppert [19] to the multivariate setting. This model will overcome the issue of inconsistent solutions arising in the FPCA approach and produce efficient results in a sparse setting. To estimate the MFAR model, we firstly form a two-level hierarchical model using the observed values. For computation, we later represent the hierarchical model as a dynamic linear model (DLM). We use observation equations to discretize the functional data and obtain measurement error. The evolution equations take into account the correlation among the functional time series and define the process model. The latent process in our case is the multivariate functional autoregressive model of order p, written as MFAR(p). A dynamic functional factor model further specifies the dynamic innovation process. We employ an efficient Gibbs sampling algorithm to perform inference and prediction for multivariate time series.
While studying the cross-correlation among the bivariate functional time series, a natural question is whether there exists a directionality of information. We are interested in analyzing the causal relations among a set of variables. As discussed by Granger [20], the first time series X is said to have a causal influence on the second time series Y if the prediction of second series Y could be improved by the inclusion of past measurements from the first time series X. Later, classical Granger causality was extended to different types of time series model. Boudjellaba, Dufour, & Roy [21] derived the necessary and sufficient condition for non-causality in stationary vector autoregressive (VAR) process. Comte & Lieberman [22] investigated second-order non-causality in the framework of the VARMA process with generalized autoregressive conditional heteroscedasticity (GARCH) type errors. This unified treatment of first and second-order causality is widely used with some modification in various fields (see, e.g., Hafner & Herwartz [23], Woźniak [24], Chen, Marron & Zhang [7]). For most of these models, restrictions on parameters are linear functions. Droumaguet & Woźniak [25] applied Markov-switching VAR models for which the restrictions for non-causality in variance are a non-linear function of parameters and used the Bayesian approach to test non-causality. Allen & Hooper [26] used Granger causality generalized measures of correlation (GcGMC) to analyze the causal relations between VIX and S&P 500.
Above mentioned literature has applied some modified versions of Granger causality to multivariate time series. Saumard [27] investigated the causality of functional time series by testing the equality of covariance operators for dependent processes. Recently, Shang, Ji & Beyaztas [28] extended the GcGMC to bivariate functional time series. This paper will address the Granger causality for functional time series utilizing the Bayes factor technique introduced by Kass & Raftery [29]. Bayes factor provides a way of evaluating evidence in favor of the null hypothesis; in our case, it tests the non-causality. Following Geweke [30], we use a sample of the posterior distribution to evaluate the model's marginal likelihood, estimate the integral over the posterior sample, and calculate the modified harmonic mean (MHM). Using the marginal density of data of the two models under consideration, we calculate and use the Bayes factor to choose the preferred model. The rest of the paper is organized as follows. In Section 2 we illustrate the MFAR model of order p and its representation in a dynamic linear model (DLM), which is further used for inference and prediction. We then generalize the concept of Granger causality to functional time series. Section 3 provides the methodology followed in this paper including a detailed explanation of the Bayes factor and the method of obtaining the same. Section 4 presents the empirical results of checking the causal relationship among the bond yields of different countries and effect of meteorological factors on air pollution level. Section 5 presents the concluding remarks.

The Model
Suppose we have K functional time series defined as {Y n t }, n = 1, · · · , K on some where Y n t ∈ L 2 (T ) for n = 1, 2, · · · , K.

Multivariate Functional Autoregressive Model
Multivariate functional autoregressive model of order p, written as MFAR(p): where is a vector of innovation functions with each ε n t being a sequence of independent mean zero random functions with E || ε n t || 2 < ∞ for n = 1, · · · , K and Ψ Ψ Ψ l = [Ψ i, j l ] K i, j=1 where each Ψ i, j l is a bounded linear operator on L 2 (T ) for i, j = 1, · · · , K and l = 1, 2, · · · , p. We assume that each function Y n t (τ n ), n = 1, · · · , K is observed at unequallyspaced points τ n ∈ T . We restrict to p = 1 and consider only integral operators for better interpretability and computational convenience. In practice, the functional observations {Y t } are observed via discrete samples of each curve with some measurement error. Suppose we observe y n i,t ∈ R sampled with measurement error v n i,t from {Y n t } ∈ L 2 (T ), n = 1, · · · , K: where τ n i,t for i = 1, · · · , m n t are observation points of Y n t and v n i,t is a mean zero measurement error with finite variance. Assume α α α t = Y Y Y t − µ µ µ, then two level hierarchical model is given as: for t = 2, · · · , T, where ψ n,m are integral operators and we assume that {v v v i,t } and {ε ε ε t } are mutually independent.
The presence of measurement error in model (3) increases the estimation error of Ψ Ψ Ψ and produces inefficient forecasts. The following proposition illustrates that MFAR(p) model for the observables is inappropriate.
and suppose we observe y y y t = Y Y Y t + v v v t are independent white noise processes. Then the observable process {y y y t } follows a multivariate functional autoregressive moving average (FARMA) process of order (p, p). For proof see Kowal et al. [19].
To overcome this model mis-specification, we separate the observed data into functional process and measurement errors in the dynamic linear model (DLM).

Dynamic Linear Models for MFAR(p)
For practical implementation of model (3), we select a finite set of evaluation points for each of functional time series, T n e = {τ n 1 , · · · , τ n M } ⊂ T for n = 1, · · · , K. We approximate the integral in (3) using quadrature methods with Q n,m for n, m = 1, · · · , K, as a known quadrature weight matrix. Let Z n t be m n t × M, n = 1, · · · , K, incidence matrix that identifies the observation points observed at time t for n-th functional times series. We can formulate the dynamic linear model (DLM) corresponding to the hierarchical model (3) as: where y y y t = [y n i,t ], is a discrete sample of K-variate functional time series {Y Y Y t } for i = 1, · · · , m n t and n = 1, · · · , K. Here, Other matrices are defined as: Model (4) can be extended to multiple lags to the MFAR(p) by replacing the evolution equation with α α α t = ∑ p l=1 ψ ψ ψ l Q Q Qα α α t−l + ε ε ε t .

Granger Causality in Functional Autoregressive Model
Here, we investigate the causal relationship between two functional time series data Y t and X t , t ∈ T observed on unequally spaced points τ ∈ T where T ⊂ R is a compact index set. Firstly, we model Y t as functional autoregressive model of order 1, written as FAR (1), and partition the linear projection of Y t on Y t−1 and X t−1 to account for the Granger causality as In second model, we model Y t as linear projection of Y t−1 only as where Y t and X t ∈ L 2 (T ) and µ Y and µ X are mean of {Y t } and {X t } respectively.
as zero mean independent random functions with covariance matrices K ε and K ε respectively. In model (M1), Y t depends on both Y t−1 and X t−1 and is considered as unrestricted model while model (M2) is considered as restricted as Y t depends only it's past values. In this setting, X is said to Granger cause Y if model (M1) is significantly better than model (M2) at forecasting future values of Y . Model (M1) and (M2) are in MFAR form, but for the computation purpose, we need to convert the MFAR models into DLM form as described in equation (4). Let Θ 1 and Θ 2 be the vector of parameters corresponding to model (M1) and (M2), respectively. Now,

Methodology
We hypothesize that the unrestricted model (M1) yields better prediction in contrast to the restricted model (M2). The null hypothesis indicates no additional contribution from the past of X to the future evolution of Y once the past of Y has been take into account. The alternative is that there is an additional contribution. In the univariate setting the usual way to test for Granger causality is to look at the forecast error sum of squares from each mode land compare them through an F test. For multivariate data, this can be extended by taking the trace or determinant of the covariance matrices under each model. As the dimension of the data becomes large, the F test becomes more unreliable as the estimation of the covariance matrix in higher dimensions is not efficient. For functional data, which is infinite dimensional, this procedure is unusable. Hence we look into Bayesian procedures which concentrate on the dimension of the parameters and not of the data. The Bayes factor is a central quantity of interest in Bayesian hypothesis testing. The ratio of the marginal data densities for both hypothesis-specific models is known as the Bayes factor. Let p(Y t | M) be the marginal data density for model M; then the Bayes factor would be Bayesian hypothesis testing begins by specifying different prior distributions on the parameters involved in the models. We also have the likelihood for the data given the prior values of the parameter. Using these, the crucial step is the computation of marginal data densities, p(Y t | M) for both models. For computing this value, we apply the Modified Harmonic Mean (MHM) method. Let Θ be a K-variate vector of parameters corresponding to the model. Firstly, obtain the sample of draws {Θ (i) } S i=1 , where S is the number of samples, using the posterior distribution of the parameters p(Θ | Y t ) for each model. Now, the marginal density of data is computed as: where, L (Y t ; Θ (i) ) is the likelihood function and p(Θ (i) ) is the prior density function of parameters obtained by substituting sample of draws {Θ (i) } S i=1 . Also, h(Θ (i) ) is a K-variate truncated Normal distribution with mean equal to posterior mean of Θ and variance set to the posterior covariance matrix of Θ.
In the rest of this section we elaborate the various steps involved in this procedure. In sections 3.1 and 3.2 we present the assumed prior distribution and the likelihood of the data under the two models. In section 3.3 we describe the method to draw samples from the posterior distribution of the parameter. Finally in section 3.4 we present the algorithm to compute the Bayes factor.

Prior Distribution
We specify the prior distributions for each of parameter in the vector Θ = (µ, Σ ν , Ψ, ε t ) below.

Likelihood Function
We use the dynamic linear model described in equation (4) to form the likelihood function for Y t . At initial time i = 1, we have Now, for i = 2, · · · , n, one-step ahead distribution of Y t is given as: Priors for each time is given as:

Posterior Sampling Distribution
We use Gibbs sampling algorithm to obtain the posterior sampling distribution. Firstly, we initialize the DLM (4) parameters and later we present the sampling distributions.
Step 3: We obtain the covariance matrix K ε ε ε and its inverse K −1 ε ε ε using the similar approach as in Step 4 of initialization.

Algorithm and Interpretation of Bayes Factor
Algorithm to compute Bayes factor: (µ, Σ ν , Ψ , ε t ) to obtain the marginal data density corresponding to model M2, p(Y t | M2). 5. Compute Bayes factor B 12 using the marginal data densities of model M1 and M2 and equation (5).

Interpretation of Bayes factor:
Kass & Raftery [29] provided the interpretation of the Bayes factor values in terms of weighing evidence against the null hypothesis, presented in Table 1. A positive logarithm of the Bayes factor is interpreted as evidence in favor of the unrestricted model. Symmetrically, a negative logarithm of the Bayes factor indicates that the data prefer the restricted model.

Simulation study
In this section, we compared the predicted results of simulated functional time series using the proposed method for MFAR relative to FAR. We are particularly interested in one-step and five-step forecasting and finding the causal   each kernel is rescaled to pre-specified squared norm, C Φ = Ψ 2 (τ, u)dτdu, with C Φ < 1. For FAR(1) we select C Ψ = 0.5 and for MFAR(1) we select C Ψ Y ,C Ψ X = (0.5, 0.1). Thus, functional time series {Y t } generated from MFAR(1) kernels is correlated to functional time series {X t } while {Y t } generated from FAR(1) kernel is independent of {X t }. We use smooth Gaussian process for innovation process ε t . In our experiment, we consider dense sampling design using m t = 30 equally-spaced observation points on [0, 1] for all t, however the model works well with sparse sampling design also. Data sets generated using FAR(1) kernel and MFAR(1) kernels are applied to proposed model (MFAR) as well as to univariate model proposed by Kowal et. al. (FAR) [19]. Table 2 shows the one-and five-step root mean squared forecasting errors (RMSFEs) for both the models. When {Y t } is dependent on {X t },we obtain lower root mean squared forecasting errors for the MFAR model than the FAR model suggesting MFAR performs better than FAR in forecasting the functional time series.
To acknowledge the causal relationship between {X t } and {Y t }, we calculate the Bayes factor as described in section 3.4. We consider three different simulated cases: (1) {Y t } is generated from MFAR and {X t } is generated from FAR, and we check whether the {X t } series causes {Y t } series, (2) {X t } is generated from MFAR and {Y t } is generated from FAR, and we check for causality in the opposite direction, and (3) {X t } and {Y t } both are generated from FAR. Table 3 presents the Bayes factor B 12 and logarithm of Bayes factor ln(B 12 ) values. Based on ln(B 12 ), we conclude that there is statistical evidence in favor of a causal relationship in first two cases and no such relationship exists when {X t } and {Y t } are independent.
We now repeat the analysis presented in Table 3 forty times for each of the dependent and independent cases and compute the Bayes factor. The boxplot of the logarithm of these Bayes factors is presented fin Figure 1. It can be seen that the distributions are well-separated. Therefore the method is useful in detecting whether the data arises from the dependent or the independent model.

Application to Real Data
We apply the proposed method to two real datasets. The first application is to daily data of yield curves over different maturities for two countries. For the second application, the response variable is still functional but the predictor variable is multivariate. The response is the daily curve of particulate matter in a city while the predictors are several meteorological variables.

Testing causality between yield curves of two countries
To assess the performance of the proposed model, we examine the causal relationship between the yield curves of two countries: USA and UK. The yield curve depicts the relationships of bond yields with different maturities. Yield curves have a significant impact on the money supply within the economy and thus give insights into the country's economic condition.

Testing whether meteorological factors cause air pollution level
Numerous pieces of literature have been forecasting air quality using meteorological features (see, e.g., Aneja et al. [31], Kumar & Goyal [32], Meng et al. [33]). Here, we investigate the causal effect of meteorological factors such as temperature, humidity, wind speed, and solar radiation on air pollution concentration. To compute particle pollution levels, we use particulate matter with diameters 10 micrometers and smaller, denoted by PM10. For our study, we used daily data from Jan 2015 to June 2020 in Delhi, India, consisting of 2008 days. However, after data cleansing we are left with 1670 days of data. We downloaded the data from official Central Pollution Control Board of India website. We take Y t as PM10 in Delhi for each hour of a day; that is, [0,23] is the support of each function. Figure 2 provides with the surface plot depicting the pollution level in Delhi over the years 2015 and 2020. For X, we take the daily average of three meteorological factors: temperature, relative humidity, and wind speed obtained from meteoblue.com. In this case, X is a vector instead of a function. Using X t and Y t , we obtain unrestricted and restricted models and apply the algorithm to compute the Bayes factor as mentioned in section 3.4. In this case, the estimated parameter space Θ is of length 361. The logarithm of the Bayes factor is obtained as ln(B 12 ) = 78.61206652. Based on the interpretation of ln(B 12 ) from table (5), we conclude that there is strong evidence that meteorological factors Granger cause air pollution.

Conclusion
The paper introduces the two-level hierarchical Gaussian process for modeling time series of multivariate functional data. We examine the Granger causality among the functional responses utilizing Bayesian techniques. We compare the marginal data densities of the restricted and unrestricted model to form a conclusion on the causal relationship among the functional time series. We implement the procedure on two real datasets. For the data on yield curves we find that the yield curve of USA is not Granger caused by that of UK. Nor is the yield curve of UK Granger caused by that of USA. In the pollution example we find that daily pollution level curves are indeed Granger caused by meteorological factors. The details of simulations and applications to actual data can be accessed through the GitHub repository Bayesian-Testing-Of-Granger-Causality-In-Functional-Time-Series.