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.
80,640 characters · 28 sections · 78 citation commands
Econometric vs. Causal Structure-Learning for Time-Series Policy Decisions: Evidence from the UK COVID-19 Policies
\email{[email removed]}
\email{[email removed]}
\affil*[1]{\orgdiv{Bayesian Artificial Intelligence research lab, MInDS research group}, \orgname{School of Electronic Engineering and Computer Science, Queen Mary University of London}, \orgaddress{\street{Mile End Road}, \city{London}, \postcode{E1 4NS}, \country{UK}}}
\abstract{Causal machine learning (ML) recovers graphical structures that inform us about potential cause-and-effect relationships. Most progress has focused on cross-sectional data with no explicit time order, whereas recovering causal structures from time series data remains the subject of ongoing research in causal ML. In addition to traditional causal ML, such as score-based and constraint-based algorithms, this study assesses econometric methods that some argue can recover causal structures from time series data. They have been used for such purposes in economics, biology and other fields. The use of these methods can be explained by the significant attention the field of econometrics has given to causality, and specifically to time series, over the years. This presents the possibility of comparing the causal discovery performance between econometric and traditional causal ML algorithms. We seek to understand if there are lessons to be incorporated into causal ML from econometrics, and provide code to translate the results of these econometric methods to the most widely-used Bayesian Network R library, bnlearn. We investigate the benefits and challenges that these algorithms present in supporting policy decision-making, using the real-world case of COVID-19 in the UK as an example. Four econometric methods are evaluated in terms of graphical structure, model dimensionality, and their ability to recover causal effects, and these results are compared with those of eleven causal ML algorithms. Amongst our main results, we see that score-based causal-ML methods often recover more identifiable effects, but typically via much denser graphs. Methodologies based on the idea of shrinkage can surface plausible intervention targets, yet no single class dominates across all metrics. Overall, the evaluation suggests complementary strengths. Econometric methods provide clear rules for temporal structures, whereas causal-ML algorithms offer broader discovery by exploring a larger space of graph structures that tends to lead to denser graphs that capture more identifiable causal relationships.}
Focusing on associations can be misleading for policy design and evaluation; we need to determine if changes in one variable cause changes in another. Traditional Machine Learning (ML) models are unsuitable for answering this question, as they focus on associations and cannot simulate the impact of hypothetical actions to inform decision-making. One must rely on influential or causal assumptions that these models do not capture.
Causal Machine Learning (Causal ML) focuses on unsupervised learning algorithms that aim to recover causal relationships from both observational and interventional data. In this study, causal ML refers specifically to causal structure learning algorithms, which focus on learning dependency structure amongst available variables, with additional causal assumptions required to interpret that structure as a causal model. This can be thought of as a map of cause-and-effect relationships of these variables. The two traditional classes of structure learning algorithms are score-based and constraint-based algorithms. Score-based learning focuses on searching the space of possible structures and optimising a score that reflects how well the learnt structure fits the data, whereas constraint-based learning uses conditional independence tests to infer causal structures. Hybrid algorithms that combine these two classes of learning are often viewed as a third class of structure learning.
Over the past few decades, there have been substantial methodological improvements in causal structure learning algorithms. However, their real-world impact remains more limited than that of traditional associational ML approaches. This is because discovering causal relationships from observational data remains extremely difficult and, unlike associational ML, small errors in causal discovery can lead to large negative repercussions in causal inference. This might explain why the literature on policy recommendations remains predominantly from the field of econometrics. That is, while other fields have undoubtedly contributed to developments in optimal decision-making, economics remains one of the most deeply engaged disciplines in these topics.
bib1 is an excellent example of a widely influential discussion on methodological approaches in econometrics that favours research designs that resemble Randomised Controlled Trials (RCTs). However, bib1 also mentions the developments in the ML field and their potential. Instrumental variables, regression discontinuity designs, and difference-in-differences methodologies remain the preferred tools for causal inference amongst econometricians. These methods have produced exciting work, such as bib2, which employs a regression discontinuity design to estimate the causal effect of immigrant legalisation on the crime rate of immigrants in Italy, concluding that legalisation reduces the crime rate of legalised immigrants. However, like Causal ML, these methods are not free of limitations. bib1 summarises issues related to weak or flawed identification strategies and problems related to this analysis's inability to generalise to other settings.
The limited development of Causal ML for time-series data presents a valuable opportunity to evaluate how well these methods support policy decisions in a real-world setting under such conditions, and to compare them with econometric techniques that are routinely used for that purpose. This study will address causality in time-series data, which adds a layer of complexity to this type of analysis. In the era of big data, massive amounts of data are continually generated over time, evolve with time and respond to external shocks, such as policy changes. The selected use case is COVID-19 in the UK, as described in detail in Section (ref). Stating this as our objective implies that we will be evaluating the following:
The development of Causal ML in the context of recovering cause-and-effect relationships from time-series data has been limited compared to econometrics. We were able to identify the following studies related to causal ML with time series data:
Although these studies are relevant, they aim to identify changes in causal structures. In contrast, we are interested in distribution shifts where the graph remains fixed but inputs, including seasonality, policy, behaviour-driven shifts, amongst others, alter observed distributions without implying new causal links.
Econometrics offers extensive time series approaches suitable for addressing this issue. A detailed explanation of how these methods can aid in recovering causal structures is provided in Section (ref). However, it will become clear that the time-series analysis methods presented work similarly to traditional causal ML methods. This is because a non-significant coefficient suggests that, after controlling for the other variables in a model, no evidence for dependence between the predictor and the target was found. This outcome is consistent with the predictor and the target being conditionally independent. The idea is similar to conventional causal discovery methods, which test for conditional independence between variables. However, regression-based inference can sometimes reflect spurious associations if important confounders are omitted. Another econometric method used involves searching for network modules, which essentially means identifying groups of interrelated nodes. This is what methods that recover causal structures do.
After learning the structures, the models can be converted into Dynamic Bayesian Networks (DBNs). A DBN is a probabilistic graphical model specifically designed to represent variables and their dependencies across different points in time. Similar to a standard Bayesian Network (BN), a DBN represents variables as nodes and their dependencies as directed arcs via a Directed Acyclic Graph (DAG). However, unlike a standard BN, where each variable is represented by a single node, a DBN represents each variable using multiple nodes, each corresponding to the same variable at different points in time. Arcs in a DBN are drawn between variables at successive time points whenever these variables are not independent of the past variable after accounting for all other relevant past variables.
In terms of parameterisation, discrete BNs typically represent dependencies using Conditional Probability Tables (CPTs), while continuous BNs employ conditional probability distributions. Under certain assumptions, BNs or DBNs can be interpreted causally, transforming the DAG into a Causal Bayesian Network (CBN). CBNs are important to this study because they enable us to model the effect of hypothetical interventions and estimate their effects, therefore enabling optimal decision-making under uncertainty.
Despite their utility, most of the applied work on BNs has been targeted to cross-sectional data (data collected at one specific point in time across multiple entities), or it has assumed that the data has been generated from multiple independent realisations of a process (repeated cross-sectional data pooled together and treated as a cross-sectional dataset). This is because BN structure learning algorithms are usually not designed to capture the evolution of a system over time, the variables that compose this system, or how their interrelations change with time.
To evaluate the usefulness of these methods for supporting policy decisions, we have selected the debated policy responses to COVID-19. The COVID-19 pandemic, caused by the SARS-CoV-2 virus, posed significant public health and policy challenges globally.
For example, vaccines are typically regarded as the primary strategy for managing the spread of viral diseases, and this expectation carried over to COVID-19 once effective vaccines were developed. Other actions related to managing the spread of this disease were viewed by many as temporary measures until a vaccine became available. This may be because traditional models for infectious disease epidemics, inspired by bib6, rely on the assumption that any member of the population can infect or be infected by any other. Many models disregard the actual interactions within a population and how diseases spread. However, it is most likely that this expectation persisted because the economic costs of implementing social distancing measures are high.
Since traditional models do not account for actual interactions within a population, and the cost of social distancing measures is high, the question arises: Which interventions are most effective in reducing the infection rate? COVID-19 will not be the last pandemic, so this crucial question is: since neither vaccines nor treatments will generally be readily available at the beginning of pandemics, can we utilise knowledge on ways to limit interactions within a population as a faster way to control pandemics? Modelling these effects by simulating hypothetical interventions in a causal inference framework bib7 may help support decision-making in future pandemics.
We use routinely collected publicly available aggregated data, which was widely used during the pandemic in the UK. For a detailed description of how the data was collected, please refer to bib8. The additional relevant data collected for this study were obtained from the same sources as those cited in bib8. These are bib9 and bib10. The dataset used in this study contains 46 columns and 866 rows, where each column corresponds to a variable and each data row to a daily outcome. The first column of the data refers to the time indicator, 'Date', which spans from January 30, 2020, to June 13, 2022. The dataset contains both continuous and categorical variables, as well as missing data values. Descriptions of the variables are provided in Table (ref).
Out of a total of 46 variables in the dataset, 40 variables contain at least one missing value, corresponding to approximately 86.96% of all variables. Overall, there are 5,667 missing values, representing 14.23% of the data values. Table (ref) offers a summary of missing data per variable. Key information includes:
To enable structure learning, we imputed the missing data. Because the data are time series, we used a state space model for imputation since it accounts for temporal dynamics when estimating unobserved values. It represents time series data as a system of state and observation equations. The state equation captures the evolution of the underlying “states" of the system over time, while the observation equation relates the states to the observed data.
The Kalman filter (bib11) is a recursive algorithm used to estimate the states in a state-space model. It is an optimal estimator when the noise is Gaussian. The Kalman filter proceeds in two steps:
When data points are missing, the Kalman filter will use only the prediction step to estimate the missing value and skip the update step for that time point. To enable the application of this methodology to categorical data, we created dummy variables, imputed them, and then converted them back. The reconversion (from dummy variables to categorical) takes the category with the highest imputed probability as the most likely category for the missing data point. It is essential to clarify that no constant was added to this model to ensure that a dummy variable represents each level of the categorical variable. The imputation was done in R using the na_kalman function of the imputeTS package (bib12). However, the last five rows still contained missing values after imputation. This might be explained by the earlier observation that most variables show higher rates of missingness towards the end of the time series. Those last five rows were omitted from the analysis.
Except for one of the econometric methods used in this study to perform structure learning, all require numeric data. Therefore, categories will be converted to numbers, so numeric data will be used for structure learning with all algorithms under consideration. However, we use the discretised dataset to later parameterise the BNs. This choice is due to the implementation restrictions in the BN software we worked with (bnlearn), which follows the common assumption that continuous nodes may only appear as children of discrete variables. Discretisation avoids these constraints and ensures that all nodes can be treated consistently during inference.
Given that most variables represent levels of intensity, such as indices measuring transportation usage or the number of hospitalised patients, we adopted the convention of discretising any continuous variables into three categories: low, medium, and high. We used the k-means clustering algorithm for this; specifically, the 2-means variant, implemented via the kmeans R function (bib13). The cluster centroids were sorted and used to determine breakpoint positions. The cut function is then used to categorise the original data into bins based on these breakpoints.
This section outlines the methodologies applied to learn the structure of the BNs and convert them to DBNs. We first describe the econometric methods used: the Least Absolute Shrinkage and Selection Operator (LASSO), the Least Angle Regression (LAR), James-Stein Shrinkage (JS), and the Statistical Inference for Modular Networks (SIMONE). We then describe the causal ML algorithms evaluated, categorising them into constraint-based, score-based, and hybrid methods. Finally, we present the common-knowledge graph constructed by bib8, which serves as one of the benchmarks for evaluating the performance of the algorithms implemented in this study.
Assume we have a DBN with a DAG \(G\) which describes a discrete-time stochastic process \(X={X_i(t); i= 1, ..,k; t = 1,..., T}\) taking values in \(R^k\) with \(k\) variables at \(t\) time points. If the model satisfies the following assumptions, it can be represented as a DBN:
The econometric methods model the dependence relationships in DBNs using a Vector Auto-Regressive (VAR) process. We say that a multivariate time series \(Y_t\) follows a VAR process of order \(p\), if:
Where:
\footnote{Please refer to bib14 for further details in this methodological discussion.}
bib15 focuses on applying the LASSO for variable selection in high-dimensional graphical models. bib16 first proposed the LASSO to address high-dimensional data analysis where the number of predictors far exceeds the number of observations. bib15 is responsible for enabling the extension of this methodology to produce graphical models. A DAG can represent the conditional dependencies between variables, and the process in bib15 can be employed to learn the structure of this graph. Focusing on a neighbourhood selection scheme using the LASSO for sparse, high-dimensional graphs, they showed that this method can recover the true graph with high probability as the sample size grows. The general idea behind using the LASSO is that real-world networks should be sparse, meaning that each variable is only directly influenced by a small number of other variables. Following the original LASSO, this is a constrained estimation procedure that tends to shrink some coefficients exactly to zero by applying an $L_1$ norm penalty to the sum of the coefficients. Only coefficients significantly different from zero define significant dependence relationships.
We implement it in R using the lars package (bib17) in the following way:
The implementation of the LASSO in the lars package (bib17) was used with the following hyperparameter:
LAR (bib18) aims to identify a suitable subset of predictors for linear regression, particularly in situations where the number of predictors exceeds the number of observations. This allows us to model interdependencies in multivariate time series, where potential models would generate many potential predictors. Therefore, the reason this is useful for producing graphical models in a multivariate time series scenario is equivalent to the usefulness of the LASSO. In essence, they address the same problem in a slightly different way.
The estimation procedure can be described in the following steps:
The relationship between this method and DAGs is no different from the previous method. Just like the LASSO, this method will return a sparse DAG. This is because both methodologies share the principle that when there are many potential causes, we should expect only a small number to have a significant influence on other variables. We also implement this method using the lars package (bib17) in R. The logic behind the implementation follows the steps described for the LASSO above. The hyperparameters need to be adjusted to type =“lar".
bib19 introduces a simple algorithm to transform correlation networks (which are undirected) into DAGs, which can imply causal relationships. The algorithm proposed approximates the full search over all possible DAGs based on the idea that in a causal network, when two variables are correlated, once the influence of all other nodes on these two nodes is removed (partial correlation), their direct link can be directed according to a given order. The direction of the influence is determined by the time order; this means the lagged variable must precede the response variable in time. In time series analysis, partial correlation can help understand direct relationships between two different time series by controlling for the influence of other series. This is particularly useful in our scenario of multivariate time series; isolating the effect of individual variables is critical. This is implemented in R as an efficient estimator of the covariance matrix, which can be obtained by allowing the estimated correlation coefficients to shrink towards zero and their estimated variance to approach its median. The package that enables this implementation in R is GeneNet bib20. The following are the hyperparameters for this method:
SIMONE (bib21) proposes a statistical framework (and its implementation in R) for inferring modules in networks. These modules or clusters tend to be characterised by dense connections among nodes within the same cluster and fewer connections between nodes in different clusters.
Using a probabilistic model, it essentially captures the likelihood of connections between nodes based on their module memberships. This approach also incorporates a form of regularisation to discourage overly complex modules. To find the most probable modules, they proposed an optimisation procedure that considers the regularisation to maximise the model's fit to the observed data. They tested their proposed methodology on both simulated and real-world data. It effectively recovered the true modules in simulated cases but was less successful in real-world data, where it could only reveal meaningful modules.
The R implementation of SIMONE has the following hyperparameters:
Above, we discussed the four different econometric methods employed to learn the structures of DBNs. Structures learnt by different methods tend to be highly inconsistent across different methods. To address this issue, we also perform model averaging on the four graphs learnt by the econometric methods to obtain an overall structure. The model-averaging output will be considered in the analysis alongside the four graphs learnt independently by each econometric method. The process employed for model averaging follows the same employed in bib22, and can be described as follows:
Structure learning is an unsupervised learning process. Its traditional version is mainly based on combinatorial optimisation, where most algorithms aim to find a discrete-valued adjacency matrix. This matrix represents the presence (and their orientation) or absence of edges.
Each of the three classes of learning described in Section (ref) is accompanied by its assumptions, merits, and limitations. There is no consensus regarding the preferred or appropriate class. Moreover, even algorithms within the same class of learning often produce widely different graphs. This phenomenon can be attributed, in part, to the fact that many algorithms operate on dissimilar assumptions about the input data and rely on distinct objective functions. Most of the considered objective functions are score-equivalent, yielding the same score for graphical structures belonging to the same Markov equivalence class. This characteristic implies that not all directed relationships can be recovered from observational data.
Algorithms that generate a unique DAG structure utilise objective functions that are not score-equivalent, or merely generate a random DAG from the highest-scoring equivalence class. This equivalence class is depicted by a Completed Partially Directed Acyclic Graph (CPDAG), which includes undirected edges (in addition to directed ones) that cannot be orientated from observational data alone. When faced with this issue in this study, we will use the cextend command from bnlearn (bib23) in R, which offers a consistent extension from CPDAGs to DAGs.
The algorithms used in this study are:
The constraint-based and hybrid algorithms require the user to input the statistical test that will be used to assess conditional independence between variables. The test selected is “mi-g-sh". According to bnlearn (bib23), this test is a shrinkage estimator based on the James-Stein estimator for mutual information. The significance level was left at the standard 0.05. The score-based and hybrid algorithms require the user to input a score that evaluates the quality of the structure learnt. The score selected is “ebic-g". According to bnlearn (bib23), "ebic-g" extends the BIC score to add a second penalty penalising dense networks.
Both the test and score selected follow the principle of penalising dense networks. This seems a fair way to compare these methods with the econometric methods used in this study, which also employ shrinkage. As with previous relevant studies (bib8,bib22), all other hyperparameters were used with their default settings. It remains unclear in the literature how to systematically compare results across different types of algorithms with varying kinds of hyperparameters. Given the extensive number of algorithms, we invite interested readers to check the default hyperparameters in bib23.
In Figure (ref) below, you will find a visual representation of the knowledge graph assumed in this study. Some variables were grouped to simplify the visual:
For a detailed explanation of how the knowledge graph was built, please refer to bib8. However, some variables were not included in the knowledge graph because we use more variables than bib8, and the logic of how to integrate them was not straightforward. Specifically, :
The knowledge graph used in this study is based on the one described in bib8, and reflects the relationship between variables relevant to the COVID-19 pandemic. It was selected for this study due to its relevance and illustrative value in demonstrating our methodology. bib8 clarifies the issues they faced in modelling feedback loops, and this paper builds on the static relationships as defined in that study and extends them to temporal relationships. Figure (ref) presents the knowledge graph used in this study, respecting the principles as defined in bib8. The starting point of this dynamic process is the pandemic, with infection and hospitalisation rates rising before government policies are implemented.
After running the econometric and causal ML algorithms using their respective libraries, we convert the resulting DAGs they produced into BN models using the bnlearn library (bib23). Since all of them produce a DAG structure, this conversion was straightforward. To extend these models into DBNs, we augment the dataset with lagged versions of the variables. The output of each algorithm is then analysed through five distinct evaluation methods. These are:
Assessing the benefits and limitations of the learnt structures to support policy decisions is more complex than assessing the structure alone. One way to validate them is to compare their predicted policy effects with the observed outcomes of those policies. Evaluating COVID-19 policies is particularly challenging due to the lack of consensus on the effect sizes, but there is a consensus on the direction of the effect. For example, it is inevitable that lockdowns reduce, rather than increase, infection rates. The learnt structures enable us to simulate various policy scenarios and determine which interventions, particularly those aimed at reducing social interactions, are most effective in lowering infection rates. If our structures consistently conclude that limiting interactions reduces virus transmission, for example, this would support their reliability for informing future policy decisions. Specifically, we will evaluate the ability of the generated structures in the following way:
The first part of this section describes how the order of the VAR process is estimated; i.e., the time-series process. Then, we analyse the normality of residuals of the VAR model and test for serial correlation. We discuss the impact of these traditional assumptions in our analysis. These diagnostic checks are essential for understanding the reliability of the estimation process and for determining what valid inferences we can make from the VAR model.
We used the VARselect function from the vars package in R (bib35) to select the optimal VAR process order \(p\) based on the following criteria:
For each model order from 1 to the maximum selected, the VARselect function fits a VAR model and calculates these information criteria. The order that minimises each criterion is then chosen. Initially, we did not set a maximum lag. However, Table (ref) shows that a maximum lag of 17 was used. We discuss this below.
Each criterion has strengths and weaknesses, and there is no correct choice for all circumstances. Results are expected to be compared across criteria, and additional diagnostic tools and domain knowledge will be used to finalise the model order selection.
AIC, HQ, and FPE find the optimal order of the VAR model at 17 lags, while SC's optimal order is one lag. We ran VARselect with a maximum possible order of more than 17 on the COVID-19 case study (meaning we ran VARselect with a maximum possible order of 18, 19, and so forth), and all of the metrics started choosing the maximum possible lag as the optimal lag (meaning all the metrics started choosing model orders of 18, 19 and so forth). This is a problem, since this parameter explosion would make the estimation of this VAR model unreliable or intractable. Figure (ref) helps us to visualise the issue. There are three possible explanations for this phenomenon:
On this basis, we proceeded with a VAR model of order one (i.e., with one lag). This would make the model more parsimonious, especially in the presence of 46 variables and 866 observations. The more lagged variables we create from the original variables, the more the number of variables will converge to the sample size. Given that generating the first lag of all our variables gives us 92 potential variables to be included in a model, we assume \(VAR(1)\) to be the optimal choice. While this might seem odd in the context of supporting policy decisions in the COVID-19 pandemic, looking at the Partial Autocorrelation Function (PACF) of New Infections shows a significant positive partial autocorrelation at lag 1. Although further lags remain significant, the decrease in the size of the PACF is considerable after lag 1, which measures the correlation between the current observation and previous lagged observations, controlling for the lags in between them. For lag 1, the PACF is equal to the Autocorrelation Function. This provides a better understanding that, although we might expect the impact of policies to take longer than one period, the last period holds most of the information we need for this analysis. This makes sense since changes will be gradual, and the number of new infections today will be closely related to the number of new infections yesterday. The PACF for New Infections is presented in Figure (ref) below.
The Jarque-Bera test (bib40) is a goodness-of-fit test that relies on the assumption that, under the null hypothesis, both the skewness and kurtosis of the population from which the sample is drawn are those of a normal distribution (i.e., skewness = 0 and excess kurtosis = 0). The test statistic \(JB\) is computed as:
Where \(n\) is the sample size, \(S\) is the sample skewness, and \(K\) is the sample kurtosis. If the sample skewness and excess kurtosis significantly deviate from zero, there is evidence to reject the null hypothesis in favor of the alternative that the data do not come from a normally distributed population. In R, this is implemented via the jarque.bera.test function of the moments package (bib41). These results are presented in Table (ref).
As shown in Table (ref), the p-values for all the series are very close to zero, leading to the rejection of the null hypothesis in all cases. This suggests that the residuals of all the variables in our VAR model do not follow a normal distribution, and this violates the assumption of the VAR process in that the error is normally distributed. However, the properties of Ordinary Least Squares (OLS) used to estimate the VAR model, on consistency and unbiasedness, do not rely on the normality assumption. This is sufficient for this study, as the existence of edges in the DAGs is derived from the coefficients of the VAR model. Since we have consistency, we know that as the sample size tends to infinity, the estimator of the coefficients converges to the true value of the parameters. Since the estimators are unbiased, we know that, on average, the estimators do not overestimate or underestimate the true values of the parameters. Therefore, we consider the coefficients obtained to be reliable. However, in the case of the JS method, failing normality implies that the Student's t-test no longer holds, and this is a limitation where we should expect JS to produce some false positive edges.
The Breusch-Godfrey test (bib42 and bib43) is used to detect autocorrelation in the residuals of a regression model. The following is a description of the process of the Breusch-Godfrey Test:
In R, this test is implemented equation-by-equation of the VAR model as follows:
The loop automates this process for each variable in the VAR model, testing each one for serial correlation. The results of this test are presented in Table (ref).
If the p-value is less than the significance level of 0.05, there is evidence of serial correlation in the residuals. These results suggest that the assumption about the behaviour of the residuals in a VAR process is violated. However, the estimate of the coefficients remains unbiased, which is adequate for the edge selection approach, which is based on coefficient magnitude. Still, this is acknowledged as a limitation.
The first subsection focuses on structure learning and compares the learnt structures to the knowledge graph, as well as metrics that describe the fit and complexity of the learnt structures. The subsequent subsection discusses policy evaluation, focusing on understanding the quality of the structures learnt concerning their application in decision-making.
The numerical results are summarised in Table (ref). The model-averaged structure (“Average") includes any edge that appears in at least two of the four algorithms. \footnote{The LASSO and LAR could not be parametrised because of the sizes of CPT that generated a bnlearn error. Therefore, only coefficients greater than 0.4 instead of 0 were considered. The model-averaged structure could only consider adding edges that appear in at least two of the structures for the same reason.} This approach increases the number of edges, and does not achieve the best BIC score amongst the econometric methods; SIMONE tops BIC as well as the LL score. \footnote{It is important to state that after using the cextend command from bnlearn (bib23) in R, GS edges could not be oriented to construct a DAG. Therefore, no results from GS are shown.}
SHD counts the number of changes needed to turn a learnt graph into the knowledge graph. These changes involve adding, removing, or reorientating an edge. The results show high SHD values for many algorithms, similar to the findings of bib8. For example, the structure learnt by LASSO contains 161 edges, yet its SHD relative to the knowledge graph is 246. This suggests that there is limited overlap between the LASSO-derived structure and the knowledge graph; essentially, to convert the LASSO graph into the common-knowledge structure, one would have to remove most of LASSO's edges and then add most edges from the knowledge graph. A strong conclusion that can be derived from these results is that the knowledge graph does not agree with the learnt graphs. However, Figure (ref) shows there are considerable differences in the graphs learnt by the econometric methods as well.
The number of free parameters varies noticeably across the algorithms, as expected given their substantially different learnt structures. Although the averaged model is more complex than SIMONE, its fit is weaker. Within the econometric set, SIMONE offers the best fit. The average model only improves over LASSO, LAR, and JS, but does not surpass SIMONE. At the lower-performing end in Table (ref) is the knowledge structure with a relatively weak LL score. This is not unexpected, as the learning algorithms explicitly optimise fit to the data, whereas the knowledge graph was not derived from data fitting.
The causal ML algorithms show SHD values that are mostly lower than those of the econometric methods, except for the two score-based algorithms HC and TABU. This is due to these score-based methods producing considerably more edges than the other algorithms, leading to substantially more complex structures than any of the other methods tested in this study, but with considerably better LL scores too. Conversely, the constraint-based methods learnt relatively sparse graphs, yet still attained competitive BIC and reasonably good LL scores.
The econometric methods are at a disadvantage in terms of BIC and LL scores, and this might be because they inherently enforce temporal constraints. In our VAR(1) specification, each node at time $t$ may only have parents at time $t-1$, so the resulting DBNs contain only past-to-present edges. This consideration becomes crucial in the next section, where we discuss the utility of these methods in informing policy decisions, showing that if an algorithm cannot discern the correct temporal order of relationships, its value for policy analysis is limited.
Simulating interventions in CBNs is equivalent to manipulating a variable in the model. In this study, this is equivalent to, for example, a policy intervention aimed at reducing one of the mobility indices and estimating the impact of this action on future cases of COVID-19. We use Pearl's bib7 do-operator to simulate interventions for causal effect estimation, where \(P(Y|do(X=x))\) represents the probability of observing \(Y\) if we manipulate \(X\) and set it to the state \(x\). We use the R package causaleffect (bib46) to compute causal effects, which implements bib47 and bib48, and implement this using bnlearn (bib23) and its cpquery function (bib23), which is based on Monte Carlo particle filters.
Since the data is discretised into three different states for continuous variables, we compute the Average Causal Effect (ACE) as follows:
The ACE represents the average change in the probability of an outcome (\(Y\)) as treatment (\(X\)) changes from a baseline state, and averages the effects across the transitions of the treatment variable. When we compute $P(Y \mid \text{do}(X = x))$, $X$ refers to the value of the intervention variable at time $t-1$ and $Y$ to the outcome at time $t$. The ACEs reported below, therefore, represent the average one-step-ahead.
Equation (ref) calculates the ACE of lowering \(X\) from a baseline set to 3, where 3 identifies the state associated with the highest values these variables can take. Here, \(X\) represents the variables related to population interactions and \(Y\) represents the variables associated with the infection rate defined in Section (ref). $X$ will be set to 3 because we want to understand if changing \(X\) from a high level of population interactions to low levels of population interactions decreases the likelihood of high levels of \(Y\).
Table (ref) below presents the results related to the evaluation metrics of the number of causal effects identified and the number of times the direction of the effect matches that of common knowledge. Table (ref) presents the average causal effects that are identified and match the direction of knowledge.
As shown in Table (ref), the only algorithms that identify a substantial number of causal effects are HC and TABU (27 effects each), with 14 and 12 of those effects respectively pointing in the direction suggested by common knowledge of the pandemic (e.g., reduced population interactions imply fewer COVID-19 cases). JS identifies 2 effects, but both match the expected direction of the effect. HC and TABU yield the most identifiable effects, explained by the denser graphs they produced, and this trade-off between identifiability and model sparsity warrants further investigation. However, unlike the econometric methods, these causal-ML algorithms were not constrained to respect the chronological order of the time series. Although the learnt graphs are acyclic in the static sense, they might display edges that are incompatible with the underlying temporal order. Therefore, Table (ref) reports Average Causal Effects (ACEs) only for those effects that satisfy all of the following: (i) they are identifiable in the learnt graph, (ii) their direction is consistent with prior/common knowledge, and (iii) they are obtained from a DBN derived from an econometric method that enforces past-to-present temporal structure (here, the JS-based model).
Under these criteria, two such interventions appear in our dataset: lowering the Citymapper journeys index and lowering OpenTable restaurant activity, each of which reduces the probability of high reinfections. Interestingly, TABU and HC also identified these effects. The effects are similar in size to those of JS. For Citymapper journeys, they stand at $-0.02$ for TABU and $-0.01$ for HC. For OpenTable restaurant activity, they stand at $-0.003$ for TABU and $-0.001$ for HC. We view this as a robustness check.
These effects are modest in magnitude but align with the expected mechanism, which suggests that reduced out-of-home mobility and indoor social activity lower transmission risk. However, there appears to be consistency across algorithms regarding the effects JS identified. This suggests that amongst the COVID-19 interventions evaluated, limiting these interactions are the most effective at reducing infection rates in the UK data. This conclusion is not easy to unravel. A crucial question here, for example, is why they are undertaking these journeys. There is evidence that although households accounted for 6% of the contacts, they were responsible for 40% of COVID-19 transmissions (bib49). Additionally, individuals at greater distances require longer exposures to contract contagion (bib49). This, paired with the OpenTable restaurant results, suggests that travel involving close contact with familiar people increases the risk of contagion, since such interactions are more likely to occur in close proximity.
This study examines how causal relationships, rather than mere associations, can be identified or assumed and used for policy evaluation. It compares econometric time-series methods with causal ML approaches to assess their effectiveness in causal discovery and causal modelling. While our results show that econometric techniques offer strengths that help address the current limitations of Causal ML in time-series contexts, both approaches present trade-offs, and neither can yet be considered universally superior.
Using an evaluation strategy that includes metrics related to learning the causal structure as well as the success of these structures in computing the direction (i.e., positive or negative) of the causal effect to support policy decisions, we found that:
Across methods, the clearest and most consistent policy signal concerned reducing journeys, thereby lowering the risk of reinfection. Although effect sizes from one day to the next were modest, this aligns with the mechanism that sustained, close-contact interactions drive transmission. For decision-makers, models that encode temporal precedence and penalise complexity seem more reliable for stress-testing interventions.
Limitations include non-normal errors, pockets of residual autocorrelation, and potential omitted variables point to model misspecification risks. Discretisation (required for DBN parametrisation) may attenuate the signal, and the lack of universally imposed temporal constraints complicates head-to-head comparisons. The knowledge graph was not designed for dynamic settings, further limiting SHD interpretability.
This study's comparison of methods produces results that could guide further work. The econometric methods failed to deal with the complexity of this dataset and the omitted variables. It is also clear that more work is needed to develop methodologies that can more accurately identify cause-and-effect relationships to inform policy decisions in time-series settings. This is made clear by the fact that while it is common knowledge that preventing human contacts prevents pandemics from spreading, many of the algorithms tested failed to detect these relationships.
\backmatter
The data used in this study have been submitted alongside the paper. Code, replication materials, and additional documentation are available in the accompanying GitHub repository:
\url{https://github.com/br1pa/econometric-vs-causal-time-series}
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
This article does not contain any studies with human participants or animals performed by any of the authors.
Bruno Petrungaro: Conceptualisation, Methodology, Software, Analysis, Writing-original draft preparation, review and editing. Anthony C. Constantinou: Writing - Review & Editing, Supervision.