Copula Models of COVID-19 Mortality in Minnesota and Wisconsin

By inergency On Nov 16, 2023

[ad_1]

1. Introduction

COVID-19 was first detected in late 2019, and by early 2020, infections were recorded in most countries around the world, leading to one of the largest pandemics in a hundred years. The pandemic initially manifested as excess mortality in many countries and regions, and was not necessarily attributed to known diseases. When the cause and character of the pandemic were revealed, it became easier to attribute known mortality cases as COVID-19-caused or COVID-19-related. Both overall and excess mortality on one hand and COVID-19-attributed mortality on the other hand should be analyzed in order to observe the course of the pandemic and predict its future development (Du et al. 2020).

The primary objective of this study is to develop the framework which will make it possible to establish possible association of pandemic mortality between geographically close regions and use it for short-term mortality forecasts. Our approach is to analyze adjacent states of the U.S. using publicly open data provided by state-level government and public health offices. We will illustrate the suggested framework with the joint COVID-19 mortality model for two upper-midwestern states: Minnesota and Wisconsin. The fully parametric copula model can be used for short-term forecasting and simulation of the mortality pattern development.

Mortality data can be obtained in two ways: the first way, which we use in this paper, is to analyze COVID-19-attributed mortality using open-source data provided by state-level health organizations. This leaves open the issue of proper attribution, and will probably miss the point of indirect effects of the pandemic (deteriorating or overwhelmed health services leading to increased mortality due to cancer and heart diseases, as well as mortality caused by negative effects on the social environment, such as stress and suicides). That is why an alternative approach would be to address the overall mortality for the pandemic period or the excess mortality calculated as additional deaths over the baseline mortality measured as the long-term average (Britt et al. 2023).

The reason we chose the first approach in this paper is that along with the easy availability of weekly data, our chosen approach tackles the drawbacks of the excess mortality approach, such as the ambiguity of average calculation and the inability to sort out the causes of death. We are likely to see a higher number of suicides, and their relationship to COVID-19 is not entirely clear. On the other hand, a lower number of traffic accidents and deaths caused by competing viral infections (like the common flu) can be associated with the COVID-19-induced preventive measures (social distancing, masking, etc.).

Following the first approach of restricting ourselves to directly attributed COVID-19-related mortality, we will study two distinct but synchronized datasets for Minnesota (MN) and Wisconsin (WI) containing the weekly number of deaths reported by the states as COVID-19-related. Our approach consists of a two-stage modeling process: first, we apply time series analysis in the Box–Jenkins or ARIMA framework to take into account the week-to-week autocorrelations in Minnesota and Wisconsin data. Then, we analyze the residuals of the time series analysis (controlled for autocorrelations). What is commonly encountered in financial time series, and may be typical for similar applications, is the violation of residual normality assumption. In order to develop a fully parametric model for the joint distribution of residuals, we use non-standard Student t and skewed t models to estimate the marginals. After that, we look for an association between the neighboring states, assuming possible non-linear tail dependence, and model it via copulas. We consider both elliptic copulas (Gaussian and t-copulas) and copulas from the Archimedean family (Clayton, Gumbel–Hougaard, etc.). The emphasis is placed on appropriate copula model selection and the construction of predictive models. Comparing the data on COVID-19-related deaths to excess mortality in future might help us to discover realistic sources of mortality dependence for geographically close areas.

There are multiple approaches to COVID-19 modeling. For a review of early COVID-19 models, see (Adiga et al. 2020). Another review of epidemic forecasting was presented by the Delphi group in (Rosenfeld and Tibshirani 2021). Various ensemble models were suggested by the Center for disease control (CDC) during 2020–2023 (Johansson et al. 2020), see also (Imperial College COVID Response Team 2022). A Bayesian approach to the reconstruction of the pandemic in the U.S. was suggested in (Chitwood et al. 2022). Problems with traditional SIR models under uncertainty and model identifiability were summarized in (Piazzola et al. 2021). Some of these problems, especially in the case of complex multivariable models, include issues with data collection, overreporting, and underreporting, and cause misrepresentation.

The autoregressive integrated with moving average methods (ARIMA) time series approach to mortality modeling was introduced in (Lee and Carter 1992). A more advanced time series technique, threshold ARIMA or self-exciting threshold autoregressive model (SETAR), was used in (Watier and Richardson 1995) to analyze epidemiological time series. In the COVID-19 context, several recent publications suggested the use of time series analysis. See, for example, (Doornik et al. 2022), (Bhattacharyya et al. 2022) and (Oliveira and Moral 2021); all of these studies applied daily time series to short-term forecasting. (Sahai et al. 2020) provided daily forecasts to five large geographically separated countries. The problems associated with this approach are discussed in (Bilinski and Salomon 2021). More sophisticated prediction techniques of mortality forecasts include autoregressive neural networks (Conde-Gutiérrez et al. 2021) and model-free Euler methods (Antulov-Fantulin and Böttcher 2022). ARIMA, along with vector autoregression, was applied to COVID-19-related excess mortality data in (Britt et al. 2023). Instead of addressing covariates of COVID-19 mortality, we will concentrate on available mortality data and restrict ourselves to basic ARIMA analysis for separate geographical regions followed by copula analysis of interregional association.

Copula models have also been applied to analyze COVID-19 data, e.g., the dependence of mortality on temperature (Alanazi 2021) or other external factors (Novianti et al. 2021). A comprehensive analysis of COVID-19 data with vine copulas is offered by (Ansell and Dalla Valle 2023). A study of interstate dependence in (Kim 2022) analyzes daily mortality data collected from fifty states in the U.S.A and D.C. (Kim 2022) uses the Gaussian copula marginal regression approach (GCMR), utilizing the principal component method and ARMA structure of copula dynamic conditional correlation to model and compare local COVID-19 mortality dynamics, concentrating on the states of California, Florida, New York, and Texas.

We follow the line of research suggested by (Kim 2022) but take a different approach. First, we use weekly data to minimize the effect of seasonality and reporting delay issues. We use the deaths officially attributed to COVID-19 to avoid misinterpretation of the causes of death. Second, we apply copula analysis to the residuals of time series models, which allows us to separate cross-correlation between states from the autocorrelation effects. This approach also makes it possible to directly address non-normality (skewness and fat tails) of the residuals allowing for marginal model selection and to use essentially non-Gaussian models such as Archimedean copulas. This approach has been suggested by (Ane et al. 2008) and applied to financial applications in (Shemyakin and Kniazev 2017).

Section 2 contains a basic description of Box–Jenkins methodology as applied in the present study, concentrating on model identification and diagnostic tools. Section 3 addresses marginal distribution analysis of the residuals, Section 4 introduces copulas, and Section 5 discusses the methods of copula model selection for the joint distribution of residuals.

Section 6 contains the results of applying ARIMA with copula analysis of residuals to COVID-19-related mortality in Minnesota and Wisconsin for 2020–2023. The principal conclusions are summarized in Section 7.

2. Arima

ARIMA (Box–Jenkins Models)

ARIMA (autoregressive integrated moving averages) methodology, also known as Box–Jenkins models (Box and Jenkins 1970), also (Box 2015) stands for autoregressive models with lag (order)

p = 1, 2, \dots

for the single variable time series

Y_{t}, t = 1, \dots, T

$Z_{t} = a_{0} + a_{1} Z_{t - 1} + \dots a_{p} Z_{t - p} + ε_{t}, ε_{t} \sim N (0, σ^{2}),$

(1)

integrated with moving average models with lag $q = 1, 2, \dots$

$Z_{t} = ε_{t} + b_{1} ε_{t - 1} + \dots b_{q} ε_{t - p}, ε_{t} \sim N (0, σ^{2}),$

(2)

applied to differences $Z_{t}$ order $d = 0, 1, 2, \dots$ in the initial time series $Y_{t}$ ,

$Z_{t} = D^{d} Y_{t}, D Y_{t} = Y_{t} - Y_{t - 1}, D^{d} Y_{t} = D D^{d - 1} Y_{t} .$

(3)

Differencing with

d > 0

makes it possible to include non-stationary models in the ARIMA

(p, d, q)

framework along with stationary ARIMA

(p, 0, q)

, also known as ARMA

(p, q)

. The general ARIMA

(p, d, q)

model can be represented (Hamilton 1994) with the help of backshift (lag) operator

L Y_{t} = Y_{t - 1}

, such that

D = 1 - L

, as

$(1 - \sum_{i = 1}^{p} a_{i} L^{i}) {(1 - L)}^{d} Y_{t} = δ + (1 - \sum_{i = 1}^{q} b_{i} L^{i}) ε_{i} .$

(4)

ARIMA model selection requires the determination of parameters $p, d, q$ and subsequent estimation of $p + q + 1$ parameters $({\hat{a}}_{0}, \dots, {\hat{a}}_{p}, {\hat{b}}_{1}, \dots, {\hat{b}}_{q})$ . This estimation is usually carried out via the method of maximum likelihood or Bayesian approach. It provides both fitted values ${\hat{Y}}_{t}$ for $t = 1, \dots, T$ and forecasts ${\hat{Y}}_{t}$ for $t > T$ , calculated recursively using residuals ${\hat{ε}}_{t} = {\hat{Y}}_{t} - Y_{t}$ .

It is common to apply unit root or more general stationarity tests to determine the differencing order d; however, other considerations can be used to suggest stationarity or non-stationarity, including the behavior of ACF (auto-regressive function) and PACF (partial auto-regressive function). Then, one can use likelihood-based information criteria such as Akaike Information Criterion (AIC) (Akaike 1974) and Bayesian or Schwarz Information Criterion (BIC) (Schwarz 1978) to determine the optimal number of lags p and q. Standard definitions of information criteria in terms of likelihood function of the data x and parameter

θ

with likelihood

L (x; θ)

are: deviance

D (θ) = - 2 l n L (x; θ)

A I C = D (\hat{θ}) + 2 k

;

B I C = D (\hat{θ}) + k log n

, where

\hat{θ}

is the value of MLE, k is the number of parameters and n is the sample size. Thus, AIC and BIC favor smaller models with lower values of p and q. Notice, however, that increasing d does not formally involve introducing additional parameters; therefore, information criteria should not be used to determine the differencing order.

If the lags p and q are chosen to be too low, residual autocorrelation will manifest itself in the results of Box–Pierce and Ljung–Box tests and graphical analysis of residual ACF and PACF. If they are too high, overparameterization will typically reveal itself in the form of large standard errors of parametric estimates.

3. Marginal Distribution Analysis of Residuals

Following the methodology discussed in detail in (Shemyakin and Kniazev 2017) and already widely used in applications to finance (Ane et al. 2008), analysis of linear and non-linear (copula) dependence of the variables of interest is applied to the residuals of ARIMA models for their marginal distributions. The first step is to check the normality of the individual residuals. Technically speaking, the absence of normality will reveal a violation of basic ARIMA assumptions. However, we can treat it as an indication that, along with serial correlation, there are other factors at play, including possible cross-correlation between the variables to be further studied in Section 4.

Standard normality testing procedures include Anderson–Darling, Kolmogorov–Smirnov and Shapiro–Wilk goodness-of-fit tests. If the hypothesis of normality is rejected, other distribution models can be considered, with Student t-distribution being a logical choice. Two main effects which could be discovered in histograms and

q - q

plots of the residuals are the increased probability of large deviations (fat tails) and/or skewness bringing about asymmetry in the left and right tails. Most popular models addressing both effects include Student t-distribution model with symmetric or asymmetric tails (centered t-distribution and skewed t-distribution, see (Fernandez and Steel 1998)).

One may suggest that if we fail to find a feasible parametric model for the marginal distribution of residuals, the other possibility is to use a non-parametric model for marginals bringing about a semi-parametric copula model. For our purpose of short-term prediction and simulation of mortality development, the fully parametric approach is more attractive due to its higher prediction power. However, misspecification of the marginals may present an issue we should be ready to address.

4. Joint Distribution Analysis via Copulas

A standard (but often insufficient) measure of association between two variables is provided by Pearson’s sample correlation. In multiple applications to economics and finance, we observe Pearson’s correlation as a measure of linear dependence missing some non-linear or tail-related association. Copula models are known to provide a good tool for addressing non-linear statistical dependence.

4.1. Overview of Copula Models

The copula approach has been very popular and successful recently. It makes it possible to estimate the joint distribution of two (or more) random variables by modeling their marginal distributions first and then analyzing their dependence structure. In this paper, we will restrict our attention to pair copulas representing the case of two variables. Possible extensions using vine copulas, hierarchical Archimedean or the hierarchical Kendall’s approach are well known and their summary can be found, for instance, in (Shemyakin and Kniazev 2017).

Let X and Y be two random variables with respective c.d.f.s $u = F (x)$ and $v = G (y)$ . Then, we represent their joint distribution $P (X \leq x, Y \leq y)$ as $C (F (x)$ , $G (y) | α) = C (u, v | α)$ , where u and v are marginal distributions, C is a copula function: a binary function with special properties, and $α$ is the association parameter measuring the strength of dependence.

The main advantage of using copulas to model joint distributions relates to the Sklar’s theorem, which asserts that not only is any copula function of u and v a valid joint c.d.f. of the vector

(u, v)

, but conversely, any joint distribution function may be represented as a copula function of the marginals. This yields an attractive practical feature, making it possible to separate marginal distributions modeling (using the parametric approach suggested in Section 3 or non-parametric estimation) from parametric modeling of their association. In the simplest cases we will consider, association parameter

α

has dimensionality of 1 or 2.

4.2. Archimedean Copulas

Archimedean pair copulas can be defined through its generator, a convex continuous strictly decreasing non-negative function $φ (t; α)$ defined for $t \in [0, 1]$ , such that $φ (1) = 0$ . The Archimedean copula is defined as $C (u, v) = φ^{[- 1]} (φ (u; α) + φ (v; α))$ , where $[- 1]$ stands for pseudo-inverse. We will consider two popular one-parametric families of pair copulas, which are often encountered in applications.

C1. Clayton’s copula

$P (X \leq x, Y \leq y) = C_{1} (u, v | α) = m a x \{{(u^{α} + v^{α} - 1)}^{1 / α}, 0\}, α \in [- 1, 0) \cup (0, \infty);$

C2. Gumbel–Hougaard’s copula

$P (X \leq x, Y \leq y) = C_{2} (u, v | α) = exp \{- {[{(- ln u)}^{α} + {(- ln v)}^{α}]}^{1 / α}\}, 1 \leq α < \infty;$

For most Archimedean families, there exist simple relationships between the association parameter

α

and Kendall’s concordance

τ

defined in (Kendall 1938). Kendall’s

τ

is widely used as a non-parametric measure of association in statistical dependence studies, as a distribution-free alternative to Pearson’s correlation. We can express association

α

for both copula models above in terms of Kendall’s concordance as

C (u, v | α) = C (u, v | α (τ))

. The relationship between

α

and

τ

for model C1 is

τ = α / (α + 2)

, and for model C2 it is

τ = 1 - 1 / α

These two families are known to effectively model dependence in the tails of the joint distribution (bivariate extremes). Clayton’s copula helps us to address the lower tail dependence (positive or negative

τ

), and the Gumbel–Hougaard copula models the upper tail dependence for positive

τ

. Along with these, we can consider dual (survival) versions of these, as described in (Wifvat et al. 2020). For instance, we might consider a dual (survival) version of Clayton’s copula constructed using marginal survival functions rather than distribution functions.

C3. Dual (survival) Clayton’s copula

$P (X > x, Y > y) = C_{3} (1 - u, 1 - v | α) = 1 - u - v + C_{1} (u, v | α) .$

Clayton’s survival copula C3, along with Gumbel–Hougaard C2 is convenient when the attention is concentrated at the right (upper) tails, and therefore may be critical in mortality studies addressing epidemic outbreaks.

Extension to the case of more than two variables for Archimedean families is non-trivial when pairwise associations are different and requires the introduction of vine copulas or hierarchical structures.

4.3. Elliptical Copulas

Elliptical copulas form another important class of copula functions. One of the main advantages of elliptical copula construction in comparison to Archimedean families is its easy extension to the case of $d > 2$ marginal distributions. However, due to their definitions through the inverses of the distribution functions, one can expect a lot more technical difficulties during the parameter estimation and simulation. The most widely used elliptical copulas are Gaussian and Student t-copula families.

C4. Gaussian copula, combined with marginal distributions

u = F (x)

and

v = G (x)

of paired data

(x, y)

, defines the joint distribution

H (x, y)

of vector

(x, y)

$P (X \leq x, Y \leq y) = C_{G} (u, v | ρ) = Φ_{2, ρ} (Φ^{- 1} (u), Φ^{- 1} (v)),$

(5)

where $Φ$ is the standard normal distribution and $Φ_{2, ρ}$ is bivariate normal with zero mean, unit variances and correlation $ρ$ .

C5. Two-parametric Student’s t-copula family (Demarta and McNeil 2005) is another popular choice

$P (X \leq x, Y \leq y) = C_{t} (u, v | ρ, ν) = T_{2, ρ, ν} (T_{ν_{1}}^{- 1} (u), T_{ν_{2}}^{- 1} (v),$

(6)

where $T_{ν}$ is a t-distribution with $ν$ degrees of freedom, and $T_{2, ρ, ν}$ is bivariate t-distribution with $ν$ degrees of freedom and correlation $ρ$ . In this case, the copula parameter is two-dimensional: $α = (ρ, ν)$ .

The relationship between $ρ$ and $τ$ for all elliptical copulas is $τ = \frac{2}{π} arcsin ρ$ .

5. Copula Model Selection

Copula model selection requires both choosing a class of copulas (Clayton, Gumbel–Hougaard, Gaussian, Student, etc.) and the estimation of copula parameters for these classes. While estimation within a given class can be optimized by using maximum likelihood (MLE) or Bayesian estimators, the non-nested models obtained for different copula families can be compared via

Comparison of AIC and BIC values for representatives of each copula family may provide a good model selection tool based on the likelihood. It is worth noting that BIC carries a higher penalty than AIC for the size of the model expressed as the number of parameters, so it tends to favor smaller models. This is especially important for smaller sample sizes.

As stated above, the copula parameters for both one-parameter Archimedean and elliptical copulas can be explicitly expressed in terms of Kendall’s tau. Thus, each class of copulas contains the best candidate $C (u, v | \hat{α})$ , where $\hat{α}$ is the MLE, corresponding to $\hat{τ} = τ (\hat{α})$ . Therefore, another standard approach will require comparing the values of $\hat{τ}$ across the copula classes with the empirical concordance $τ^{*}$ calculated directly from the data. The smaller distance $| \hat{τ} - τ^{*} |$ corresponds to the better model.

To measure the goodness of fit of the model, one can define the distance between the empirical c.d.f.

F_{n}^{*} (x)

and the model c.d.f.

F (x)

. The Kolmogorov–Smirnov statistic is defined as

$d_{K S} (F_{n}^{*} (x), F (x)) = \sup ∣ F_{n}^{*} (x) - F (x) ∣$

(7)

and can be used for goodness-of-fit testing. Another good choice is the Cramer–von Mises statistic

$d_{C v M} (F_{n}^{*} (x), F (x)) = \int_{- \infty}^{\infty} {[F_{n}^{*}}^{(x) - F (x)] 2} d F (x)$

(8)

In multidimensional situations, finding critical values for these statistics might be problematic. However, if the final goal is to compare the values of

d_{K S}

and

d_{C v M}

for the best representatives (MLEs) of copula families, this relative analysis can provide a reliable tool of model selection (Shemyakin and Kniazev 2017).

Approximate p-values for the default test statistic from the Cramer–von Mises functional, as detailed in (Genest et al. 2009), can be derived either through a parametric bootstrap or a more expedient multiplier method. Using the parametric bootstrap, one can estimate dependence parameters for most copula families. With the multiplier method, rank-based strategies apply in bivariate scenarios, while only maximum pseudo-likelihood estimation is viable for multivariate cases. The computational advantage of the multiplier method demands additional programming due to required partial derivative calculations for each copula family.

Tail dependence is the characteristic that describes the behavior of bivariate extremes occurring in the tails of the distributions corresponding to the lower-left and upper-right corners of the square $(u, v) \in [0, 1] \times [0, 1]$ . The definitions of tail dependence coefficients for copulas $C (u, v)$ with marginals $u = F (x)$ and $v = G (y)$ are as follows:

Lower tail dependence coefficient

$λ_{l} = \lim_{q \to 0^{+}} P (F (X) \leq q | G (Y) \leq q) = \lim_{q \to 0^{+}} \frac{C (q, q)}{q} .$

Upper tail dependence coefficient

$λ_{u} = \lim_{q \to 1^{-}} P (F (X) > q) | G (Y) > q) = \lim_{q \to 1^{-}} \frac{1 - 2 q + C (q, q)}{1 - q} .$

One-parametric Archimedean copulas exhibit only either lower (e.g., C1) or upper (e.g., C2, C3) tail dependence, while the t-copula family is characterized by symmetric lower and upper tails. Thus, when selecting an accurate copula model for specific data, one should consider whether the data display upper and/or lower tail dependence, and whether both tails are of equal importance for the objective of the study. In our case, upper tails corresponding to pandemic outbreaks may play a more significant role. The following closed-form expressions can be used to calculate the lower and/or upper tail dependence coefficients of the selected five copula models through their implied $τ$ values: C1: $λ_{l} = 2^{\frac{τ - 1}{2 τ}}, λ_{u} = 0$ , C2: $λ_{l} = 0, λ_{u} = 2 - 2^{1 - τ}$ , C3: $λ_{l} = 0, λ_{u} = 2^{\frac{τ - 1}{2 τ}}$ , C4: $λ_{l} = λ_{u} = 0$ .

There is a special expression for tail dependence in case of t-copulas class C5 with degrees of freedom

ν

and correlation

ρ

$λ_{l} = λ_{u} = T_{ν + 1} (- \sqrt{ν + 1} \sqrt{\frac{1 - ρ}{1 + ρ}}) .$

In order to use model tail dependence as a criterion of model selection, we may want to compare the closed form results for copula families to an empirical measure of tail dependence suggested in (Brechmann 2010): empirical tail cumulation, defined for both tails for quantile q and bivariate sample

(u_{i}, v_{i}) \in {[0, 1]}^{2}, i = 1, \dots, n

$γ_{l} (q) = \frac{1}{n} \sum_{i = 1}^{n} I_{u_{i} \leq \sqrt{q}, v_{i} \leq \sqrt{q}} - q, γ_{u} (q) = \frac{1}{n} \sum_{i = 1}^{n} I_{u_{i} > \sqrt{q}, v_{i} > \sqrt{q}} - q,$

(9)

which also makes it possible to match copula models with empirical pseudo-observations $(u, v)$ . Upper tail cumulation can be defined similarly. Unfortunately, this approach, while providing some graphical insights, does not allow us to clearly distinguish different tail dependence patterns, so it has to be used cautiously.

6. COVID-19 in Minnesota and Wisconsin

We will use an example of two U.S. upper midwestern states: Minnesota (MN) and Wisconsin (WI). They share a common border and certain demographic characteristics, with substantial rural populations along with some highly populated metro areas (Minneapolis/St. Paul in Minnesota and Milwaukee in Wisconsin). It is not uncommon for residents of one state to commute to the other for work and entertainment. Traffic across the state border is substantial, utilizing seven bridges and three major highway crossings between the states. According to (O’Neill 2016), based on the U.S. census data, over 242,000 workers commuted into the Twin Cities metro area in Minnesota, with one-third of them estimated to be coming from Wisconsin. In addition, over 97,000 commuters travel to outflow jobs (including Wisconsin) from the Twin Cities metro area.

The geographical and cultural proximity of the two observed states suggests a certain communality in the pandemic development and a possible similarity in their mortality patterns. These could be due to the global and nationwide development of the disease, but also to the direct spread of the disease across the border, since no specific restrictions were applied to interstate travel. Also, Minnesota and Wisconsin share similar climate conditions, out-of-state travel patterns, and school and college calendars which may affect the pandemic spread.

However, there were also differences between how the two states and their healthcare systems addressed COVID-19. Minnesota introduced a comprehensive program of social distancing measures as early as the summer of 2020. In the meantime, Wisconsin imposed much weaker restrictions, with fitness clubs and bars remaining open for most of the summer. For winter sports, most indoor gyms and ski areas in Wisconsin were open in 2020–2021, while Minnesota had imposed substantial restrictions. Statewide and local school and college social distancing measures were also different.

We do not aim to provide a comparative analysis of the two states’ COVID-19 mortality patterns or of the effect of health policy decisions. The focus of our study is on finding dependence patterns in two states’ COVID-19 mortality numbers. For that purpose, we first apply one-dimensional time series analysis to each state’s mortality data separately. That makes it possible to filter out the major trends in disease development and autocorrelation within the time series. On the second stage, we apply copula analysis to the time series residuals. It will facilitate capturing the finer effect of the association between two time series. The purpose of this analysis is to establish statistical patterns which would allow one to use one state’s data for another state’s mortality forecast or build joint forecast models for two states.

The choice of these two upper Midwestern states is also explained by the easy availability of COVID-19 data in a similar synchronized weekly format. In both Minnesota and Wisconsin, the first COVID-19-related deaths were recorded 19 March 2020. It would be interesting to consider mortality for different age groups, but full age data were not easily available in open-source format. Further analysis is intended to be expanded to other adjacent states and different countries in the world, but in that case, the difference in public data formats may be substantial.

We retrieved raw COVID-19 mortality data from the Department of Health of MN (Minnesota Department of Health 2023) and WI (WI Dept of Health Services 2023). To perform ARIMA analysis of mortality, we needed to determine the total number of deaths per week reported by the states as COVID-19-related. Deaths with laboratory testing are attributed to COVID-19 if: (a) A positive PCR test (confirmed case) or antigen test (probable case) for SARS-CoV-2 and COVID-19 was listed on the death certificate or (b) clinical history/autopsy findings provided evidence that the death was related to COVID-19 without an alternative cause (i.e., drowning, homicide, trauma, etc.). The observation window covers 172 consecutive weeks ranging from March 2020 to June 2023. Figure 1 shows the weekly mortality rates per 100,000 in MN and WI. Two lines on the graph are similar, but it remains to be seen if this similarity can be explained by the parallel autocorrelation patterns or there is also a distinct association between two time series going beyond that.

We applied the Zivot–Andrews structural break test (Zivot and Andrews 1992) to consider a possible change of pattern in the time series. The potential break point for MN was detected at position 43; the potential break point for WI was at position 42. Both correspond to obvious return-to-school jumps in the mortality-rate curves in fall 2020. However, we chose to use entire time series for further analysis, assuming the account for minor structural breaks was not instrumental for overall modeling, but sample size was an issue. Analysis of daily instead of weekly data would increase the sample size, but would also raise issues with seasonality and misspecification of reporting dates.

We may argue that ARIMA provides an adequate tool for modeling pandemic mortality since the period of study includes both the onset and the later stages of the disease. Even with three discernable peaks during the observation period, the autocorrelation pattern seems to be at least equally important, especially for the purpose of short-term forecast. Therefore, ARIMA, or even ARMA, if the data pass the test of stationarity, presents a viable modeling choice.

6.1. ARIMA Modeling

The next step is to check the principle assumptions of ARIMA methodology for weekly mortality in Minnesota and Wisconsin. As demonstrated by the ADF (advanced Dickey–Fuller) unit root test rejected with $p < 0.01$ and the KPSS (Kwiatkowski–Phillips–Schmidt–Shin) stationarity test not rejected with $p > 0.1$ , the entire time series may be considered stationary, which is not surprising, since the observation window includes both the onset and outset of the active pandemic. Thus, the difference order in ARIMA $(p, d, q)$ is suggested as $d = 0$ .

We have tried all combinations of AR and MA components from 0 to 4 and, based on the information criteria and significance analysis of the model coefficients, we can suggest ARIMA (3,0,0) (suggested by BIC) or ARIMA (4,0,2) to describe the behavior of the time series for MN, and ARIMA (3,0,0) for WI. The final summary for MN and WI is given in Table 1. The BIC values were the smallest of all the tested lags. After estimation in ARIMA models of choice, we performed basic residual analysis tests. The Box–Pierce and Box–Ljung tests, as well as the analysis of residual ACFs and PACFs, confirmed the absence of serial correlation in the residuals. Both tests for the selected models would not reject the null hypothesis (for p-values, see the last column of Table 1).

6.2. Marginal Distribution Analysis of Residuals

Along the lines of Section 3, we began with the normality test for the residuals of MN and WI series denoted by

x_{t}

and

y_{t}

. Violation of normality would substantially affect the further results. We can graphically observe fat tails (and possibly, skewness) in the

q - q

plots. With no residual autocorrelation in

x_{t}

and

y_{t}

, they can be treated as independent samples.

Then, we performed a standard normality check, and the results of the Anderson–Darling and the Shapiro–Wilk tests indicate that the residuals are not normally distributed. Both tests yield extremely small p-values (less than

2.2 \times 10^{- 16}

), indicating strong evidence against the null hypothesis of normality and confirming what we see in Figure 2. The next model to be logically considered is the non-central Student t-distribution with p.d.f.

$t (x; ν, μ, σ) = \frac{Γ ((ν + 1) / 2)}{σ \sqrt{ν π} Γ} {(1 + \frac{1}{ν} {(\frac{x - μ}{σ})}^{2})}^{- \frac{ν + 1}{2}}$

(10)

The numerical results of fitting the non-central Student t distribution to the ARIMA residuals, along with goodness of fit, are summarized in Table 2.

Graphical support for this hypothesis is provided by the

q - q

plots in Figure 3, providing much better fit in the tails with the exception of a few outliers.

Skewness does not seem to present a significant issue, so the generalization to the class of skewed t-distribution, which we also considered

$f (x; ν, μ, σ, α) = \{\begin{matrix} 2 t (x; ν, μ, σ) T (\frac{ν}{2}, \frac{ν}{2}, α \frac{x - μ}{σ}), & if x < μ \\ 2 t (x; ν, μ, σ) [1 - T] (\frac{ν}{2}, \frac{ν}{2}, α \frac{x - μ}{σ}) \end{matrix}, if x \geq μ$

(11)

(compare (Fernandez and Steel 1998)), does not seem to be necessary. The choice of non-skewed non-standard t-distribution is confirmed as viable by the values of information criteria in Table 3 and histograms in Figure 4. This seems to present a difference from many financial time series in which substantial skewness of residuals calls for the use of a skewed t model (compare (Shemyakin and Kniazev 2017)).

In view of the test results, information criteria and graphical evidence, we will consider marginal models for the ARIMA residuals for MN

x_{t}

and WI

y_{t}

, as suggested by two respective lines of Table 2 as non-central t-distributions with c.d.f.s

F (x)

and

G (y)

6.3. Joint Distribution of Residuals

Before performing the joint distribution analysis, we can visualize the relationship between the two sets of residuals to get an idea of their dependency structure. The graphs in Figure 5 represent standardized residuals of the ARIMA model

(x_{t}, y_{t})

in the left panel and transformed values

(u_{t} = F (x_{t}), v_{t} = G (y_{t}))

using parametric estimates from Table 2 in the right panel. The left panel shows some outliers caused by both the bivariate dependence and the behavior of the marginal distributions, while the right panel illustrates lower and upper tail dependence free of the marginal effect. Both suggest the use of a copula model. Clustering of points both in the lower left and in the upper right corners may serve as a clear indicator of tail dependence.

Since both lower and upper tails seem to be present in Figure 5, we performed copula analysis for all five families described in Section 4, using the tailless Gaussian C4 as a benchmark. The best candidates in each family were selected via MLE. Some technical difficulties were encountered in the case of Student t-copula, and numerical optimization was performed using Nelder–Mead and Broyden–Fletcher–Goldfarb–Shanno methods. The data for five families are summarized in Table 4, providing information criteria values, model induced Kendall’s

τ

, calculated according to Section 4, Cramer–von Mises test p-values and tail dependence coefficients (lower in case of C1, upper in case of C2 and C3, symmetric lower and upper for C4 and C5). We can observe that the empirical Kendall’s concordance value is

\hat{τ} = 0.315

. The best fit is indicated by the bold font.

To illustrate the results of copula modeling graphically, Figure 6 depicts transformed residual data from the right panel of Figure 5 overlaid with the contour plots of fitted copula models (C1 upper left, C2 lower left, C3 upper right and C5 lower right). It appears that the clusters of points in both left lower and right upper corners of each panel (as well as some outliers in both the upper-left and lower-right corners) are better explained by the tail structure of the t-copula.

It appears that the Student t-copula is the clear overall winner. The MLE estimates of parameters for the selected model (with standard error in parentheses) are

ρ = 0.455 (0.11), η = 1.365

. As an alternative, we can consider the method of moment estimates based on empirical value of

τ

ρ = 0.475 (0.08), η = 1.365

. The main limitation of t-copulas is the symmetry of lower and upper tails. If we extend the model to skew t-copulas, which would make it possible to address two asymmetric tails, the computational and numerical challenges of t-copula modeling discussed in (Hintz et al. 2022) become even more serious. Gumbel–Hougaard and dual (survival) Clayton models also deserve attention, since they adequately address the upper tail dependence being of primary concern.

In order to provide some intuition for the tail dependence figures, Figure 7 contains graphs of empirical tail cumulation compared with the tail behavior expected of selected models. We generated 10,000 samples from each candidate copula model using a bootstrap sampling method and chose quartile intervals for the lower and upper tails (

[0.00, 0.25] & [0.75, 1.00]

The left panel illustrates the lower tail cumulation and the right panel shows the upper tails. Due to a relatively small sample size, the empirical tail cumulation curve is rather glitchy and reveals the effect of outliers. The copula model curves are shown with a 95% confidence interval for the values of model (theoretical) tail cumulation. Clayton and t-copula provide the best model for the lower tail, while dual Clayton, Gumbel–Hougaard and t-copula provide a better fit for the upper tail.

6.4. Construction of Prediction Intervals

If we consider the problem of short-term forecasting of COVID-19 mortality (with standards being one week or four weeks ahead), the point forecasts for both Minnesota (MN) and Wisconsin (WI) can be obtained by ARIMA for time horizon

{T + 1, \dots, T + 4}

using coefficients

α

from Table 1 and available mortality data

z_{1}, \dots, z_{T}

recursively as

$\begin{matrix} z_{M N, T + 1} = α_{M N, 1} z_{M N, T} + α_{M N, 2} z_{M N, T - 1} + α_{M N, 3} z_{M N, T - 2} \\ z_{W I, T + 1} = α_{W I, 1} z_{W I, T} + α_{W I, 2} z_{W I, T - 1} + α_{W I, 3} z_{W I, T - 2} \\ \dots \\ z_{M N, T + 4} = α_{M N, 1} z_{M N, T + 3} + α_{M N, 2} z_{M N, T + 2} + α_{M N, 3} z_{M N, T + 1} \\ z_{W I, T + 4} = α_{W I, 1} z_{W I, T + 3} + α_{W I, 2} z_{W I, T + 2} + α_{W I, 3} z_{W I, T + 1} . \end{matrix}$

However, in order to build prediction intervals (say, with

(1 - q) %

prediction probability), one needs to simulate a large number N of independent identically distributed pairs

(x_{T + j}^{i}, y_{T + j}^{i}), i = 1, \dots, N, j = 1, \dots, 4

, where the joint distribution of each pair is modeled by t-copula C5 with parameters

ρ = 0.455, η = 1.365

and marginals from Table 2. Then, we generate N scenarios,

i = 1, \dots, N

$\begin{matrix} z_{M N, T + 1}^{i} = α_{M N, 1} z_{M N, T} + α_{M N, 2} z_{M N, T - 1} + α_{M N, 3} z_{M N, T - 2} + x_{T + 1}^{i} \\ z_{W I, T + 1}^{i} = α_{W I, 1} z_{W I, T} + α_{W I, 2} z_{W I, T - 1} + α_{W I, 3} z_{W I, T - 2} + y_{T + 1}^{i} \\ \dots \\ z_{M N, T + 4}^{i} = α_{M N, 1} z_{M N, T + 3}^{i} + α_{M N, 2} z_{M N, T + 2}^{i} + α_{M N, 3} z_{M N, T + 1}^{i} + x_{T + 4}^{i} \\ z_{W I, T + 4}^{i} = α_{W I, 1} z_{W I, T + 3}^{i} + α_{W I, 2} z_{W I, T + 2}^{i} + α_{W I, 3} z_{W I, T + 1}^{i} + y_{T + 4}^{i} . \end{matrix}$

Now, the intervals

$(Z_{M N, T + j}^{l o w e r}, Z_{M N, T + j}^{u p p e r}), (Z_{W I, T + j}^{l o w e r}, Z_{W I, T + j}^{u p p e r}), j = 1, \dots, 4,$

with their bounds corresponding to $q / 2 %$ (lower) and $(1 - q) / 2 %$ (upper) sample percentiles of $Z_{M N}^{i}$ and $Z_{W I}^{i}, i = 1, \dots, N$ fulfil the purpose of approximate $(1 - q) %$ prediction.

7. Conclusions

We suggest ARIMA modeling as a simple tool for addressing mortality forecast for a given geographical or administrative region at all pandemic stages using only open access mortality data. The study of COVID-19-related mortality for Minnesota and Wisconsin reveals somewhat similar autocorrelation patterns, reflecting significant temporal dependence of weekly mortality on three previous weeks’ data. That should open up a possibility for a short-term point forecast of pandemic mortality based on available past data.

However, in order to evaluate errors of such point forecasting and provide at least some prediction intervals, one needs to perform a further residual analysis. This is especially important when we consider that residual distributions fail tests for normality and exhibit heavier tails. This would bring about underestimation of errors if we use techniques based on the standard normality assumption.

One more level of complexity is added by analyzing data from adjacent regions. The dependence pattern of residuals obtained for Minnesota and Wisconsin leads us to believe that additional information on weekly mortality can be provided by studying the joint distribution of residuals for several distinct areas. Copula models provide a powerful tool for studying both linear and non-linear dependence of mortality in geographically close regions. Therefore, the quality of forecast can be substantially improved by analyzing the data obtained by the neighbors. The copula model selection plays an important role. Based on our example, we can recommend the Student t-copula addressing both tails of the joint distribution and Gumbel–Hougaard and dual Clayton copulas stressing the role of the upper tail (high mortality outbreaks).

Copula models provide an opportunity to simulate joint distribution of mortality across wider geographical areas using the public data obtained for every region separately. Such simulation can generate multiple future scenarios, making it possible to estimate the probabilities of catastrophic development of a pandemic. The paper contains a description of the general methodology, as well as a specific example of COVID-19 development for two U.S. states of Minnesota and Wisconsin.

[ad_2]