key: cord-0618400-w8nnps3j
authors: Lima, Alexandre Barbosa de
title: Covid 19 and A Wavelet Analysis of the Total Deaths per Month in Brazil since 2015
date: 2020-09-01
journal: nan
DOI: nan
sha: f63686fa54043be7f951e34bb56b8a0c1b157e1b
doc_id: 618400
cord_uid: w8nnps3j

We investigate the historical series of the total number of deaths per month in Brazil since 2015 using the wavelet transform, in order to assess whether the COVID-19 pandemic caused any change point in that series. Our wavelet analysis shows that the series has a change point in the variance. However, it occurred long before the pandemic began.

The remainder of the paper is organized as follows. Section 2 presents an overview of the Continuous and Discrete Wavelet Transforms for the reader who is not familiar with the subject. If the reader is familiar with the theory of wavelets, he or she may proceed to Section 3, which presents the experimental results. Finally, section 4 presents our conclusions and and highlights some topics for further work.

The Fourier Transform (FT) of a signal x(t), t ∈ R (t denotes time), if exists, is defined as (1) in which ν denotes the frequency in cycles/second [Hz] .

Gabor [12] has shown that it is possible to represent the local spectral content of a signal x(t) around an instant of time τ by the Windowed Fourier Transform (WFT).

in which g T (t) is a window of finite duration support T and ν denotes frequency.

The WFT is a two-dimensional representation defined on the time-frequency domain (or plane) as it depends on the ν and τ parameters. The WFT would be equivalent to a kind of continuous "sheet music" description of x(t).

According to the Heisenberg's uncertainty principle [13, p.52 ], a signal whose energy content is quite well localized in time has this energy quite spread out in the frequency domain. As the window of (2) has a fixed size T , we may conclude that the WFT is not good to analyze (or identify) behaviors of x(t) occurring in time intervals much smaller or much larger than T , as, for example, transient phenomena of duration ∆t << T or cycles that exist in periods larger than T .

A wavelet ψ 0 (t) (sometimes also called mother wavelet), t ∈ R, is a function that satisfies three conditions [4] , [14] .

1. Its Fourier transform Ψ(ν), −∞ < ν < ∞, is such that exists a finite constant C ψ that obeys the admissibility condition

2. The integral of ψ 0 (t) is null:

3. Its energy is unitary: Figure 1 shows four examples of wavelet functions: Haar, Daubechies, Coiflet and Symmlet. As the name suggests, a wavelet is a 'small wave'. A small wave grows and decays in a limited time period. On the other hand, an example of a 'big wave' is the cosine function cos t, which is 'eternal', i. e., keeps oscilating up and down for all t.

The wavelet transform is a relatively new tool for the analysis of signals, given that his mathematical theory was formalized in the 1980s [15] . The wavelet transform has been originally developed as an analysis and synthesis tool of continuous time energy signals [16] , [17] , [18] , [19] , [15] , [20] .

An energy signal x(t), obeys the constraint

i. e., x(t) that obeys the constraint (6) belongs to the squared summable functions space L 2 (R). Presently, the wavelet transform has also been used as an analysis tool of discrete time signals.

There are continuous time and discrete time wavelet decompositions designate by Continuous Wavelet Transform (CWT) and Discrete Wavelet Transform (DWT).

The CWT of a signal x(t) consists of a set C = {W ψ (s, τ ), s ∈ R + , τ ∈ R}, in which

• τ is the time localization parameter,

• s represents scale, and

• ψ denotes a wavelet function, of wavelet coefficients on the continuous time-scale plane (also known as time-frequency plane) given by

denotes a dilated and shifted version of the "mother" wavelet ψ 0 (t).

The factor 1/ √ s in (7) provides all functions of the class

have the same energy (norm).

The basic idea of the CWT defined by (7) is to correlate 1 a signal x(t) with shifted (by τ ) and dilated (by s) versions of a mother wavelet (that has a pass-band spectrum). The CWT is a two parameters function. So, it is a redundant transform, because it consists on mapping an one-dimension signal on the time-scale plane.

Differently from the WFT, where the reconstruction is made from the same family of functions as that used in the analysis, in the CWT the synthesis is made with functionsψ s,τ that have to satisfỹ

So, x(t) is completely recovered by the Inverse Continuous Wavelet Transform (ICWT):

The fundamental difference between the CWT and the WFT consists of the fact that the functions ψ s,τ undergo dilations and compressions [13] . The analysis on refined scales of time (small values of s) requires "fast" ψ s,τ functions, i. e., of a small support, while the analysis on aggregate scales of time (large values of s) requires "slower" ψ s,τ functions, i. e., of a wider support. As already mentioned, the internal product defined by (7) is a measure of similarity between the wavelet ψ t−τ s and the signal x(t) on a certain instant of time τ and on a determined scale s. For a fixed τ , large values of s correspond to a low-frequency analysis, while small values of s are associated to a high-frequency analysis. Therefore, the wavelet transform has a variable time resolution (i. e., the capacity of analyzing a signal from close -"zoom in" -or from far -"zoom out"), being adequate to analyze phenomena that occur in different time scales. Note that the discontinuity produces large coefficients in its respective cone of influence, which converge to the location of the singularity.

There are two kinds of DWT:

• the DWT for discrete time signals; and • the DWT for continuous time signals.

The DWT may be formulated for discrete time signals (as it is done, for example, by Percival and Walden [4] ) without establishing any explicit connection with the CWT. On the other hand, we should not understand the term "discrete" of the DWT for continuous time signals as meaning that this transform is defined over a discrete time signal. But only that the coefficients produced by this transform belong to a subset D = {w j,k = W ψ (2 j , 2 j k), j ∈ Z, k ∈ Z} of the set C [14, p.105], [21] .

In fact, the DWT coefficients for continuous time signals can also be directly obtained by means of the integral

in which the indices j and k are called scale and localization, respectively, does not involve any discrete time signal, but the continuous time signal x(t).

Equation (11) shows that the continuous time DWT corresponds to a critically sampled version of the CWT defined by (7) in the dyadic scales s = 2 j , j = . . . , −1, 0, 1, 2, . . ., in which the instants of time in the dyadic scale s = 2 j are separated by multiples of 2 j . The function ψ 0 of (11) must be defined from a Multiresolution Analysis (MRA) of the signal x(t) [4] , [20] , [22] . Observe that the continuous time MRA theory is similar to that of discrete time.

In this paper, we decided, for mere convenience, to present the continuous time MRA version based on the spectral analysis of a "fictitious" signal {x t , t ∈ R} that is associated to the discrete time series {x n , n ∈ Z} [21] . 

The subspace V j is known as the approximation space associated to the time scale s j = 2 j (assuming that V 0 is the approximation space with unit scale).

If the x(t) projection on V j is represented by the scale coefficients

then properties 1 and 3 assure that lim

Property 4 implies that the subspace V j is a scaled version of subspace V 0 (multiresolution).

The orthonormal basis mentioned in property 5 is obtained by time shifts of the low-pass function φ j .

Consider the successive approximations sequence (also known in the literature as wavelet smooths [4] ) of x(t)

This fact illustrates the MRA's fundamental idea, that consists in examining the loss of information when one goes from S j (t) to S j+1 (t):

∆x j+1 (t) (called detail of x j (t)) belongs to the subspace W j+1 , named detail space [4] that is associated to the fluctuations (or variations) of the signal in the more refined time scale s j = 2 j and that corresponds to the orthogonal

The MRA shows that the detail signals ∆x j+1 (t) = D j+1 (t) may be directly obtained by successive projections of the original signal x(t) over wavelet subspaces W j .

Besides, the MRA theory shows that exists a function ψ 0 (t), called "mother wavelet" , that is obtained

The internal product ψ j+1,k (t), x(t) = w j+1,k denotes the wavelet coefficient associated to scale j + 1 and discrete time k and {ψ j+1,k (t)} is a family of wavelet functions that generates the subspace W j+1 , orthogonal to subspace V j+1 (W j+1 ⊥V j+1 ), i. e., ψ j+1,n , φ j+1,p = 0 , ∀n, p.

Therefore, the detail signal D j+1 (t) belongs to the complementary subspace W j+1 of V j , because

That is, V j is given by the direct addition of V j+1 and W j+1 , and this means that any element in V j may be determined from the addition of two orthogonal elements belonging to V j+1 and W j+1 . Iterating (17), we have

Eq. (18) says that the approximation S j (t) is given by

The MRA of a continuous time signal x(t) is initiated by determining the coefficients 3 u 0 (k) = φ 0,k (t), x(t) , in which k = 0, 1, . . . , N − 1, that are associated to the projection of x(t) on the approximation subspace V 0 .

Following, the sequence {u 0 (k)} is decomposed by filtering and sub-sampling by a factor of 2 (downsampling) in two sequences: {u 1 (k)} and {w 1 (k)}, each one with N/2 points. This filtering and sub-sampling process is repeated several times, producing the sequences

2 Besides, Wj+1 is contained in the subspace Vj. 3 The sequence u0(k) is obtained sampling the filter's output whose impulse response is φ * (−t) (matched filter with a function φ0(t) = φ(t)) at instants k = 0, 1, 2, . . ., i. e., u0(k) = x(t) φ * (−t) for k = 0, 1, 2, . . ., in which denotes convolution. and

The literature calls the set of coefficients [21] , [23] 

as the DWT of the x(t) signal. The reconstruction of x(t) is implemented by filtering and oversampling by a factor of 2 (upsampling) of the sequences (20) and (21), obtaining an approximation of x(t) in the subspace V 0

or

Eq. (24) defines the Inverse Discrete Wavelet Transform (IDWT).

We say that the function φ 0 (t) = φ(t) determines a MRA of x(t) according to (23) , if it obeys the following conditions:

1. intra-scale orthonormality (property 5)

in which δ m,n is the Kronecker's delta (δ m,n = 1 if m = n, δ m,n = 0 for m = n). Eq. (25) imposes an orthonormality condition at scale j = 0.

as several φ(t − k) fit in φ( t 2 ) (is a consequence of property (1) of the MRA).

Equation 27 may be rewritten as

known as Dilation Equation, n ∈ Z.

Equations 27 and 28 may be written, respectively, in the frequency domain as

and

in which Φ(ν) is the FT of φ(t) and G(ν) = n g n e −j2πνn , known as scale filter (low-pass), represents a periodic filter in ν.

As the subspace W j+1 is orthogonal to V j+1 and is in V j , we have

or

that is the Wavelet Equation.

Applying the FT to (31) and (32) we get, respectively,

and

in which H(ν) is the wavelet filter (high-pass).

Rewriting (16) in terms of the frequency domain and using (29) and (33) results the orthogonality condition

that the filter H has to obey so the family {ψ 1,k (t)} is orthogonal to the family {φ 1,k (t)}.

We may show that the condition [ 

in which L denotes the length of a Finite Impulse Response (FIR) filter g n , H(z) and G(z) denote the z-transform of sequences h n and g n , respectively, is sufficient to (35) to hold.

We say that g n e h n are Quadrature Mirrored Filters (QMF) when they are related by (36). Figure 6 shows the QMF filters frequency response related to a Daubechies wavelet of order 10 (db10). According to (28), the MRA departs from a definition (from several possible) of the scale function φ(t), that is related to the scale filter g n by (27). Eq. (36) says that the choice of a FIR-type filter {g n } implies a {h n } that is also FIR. At last, the wavelet function is determined by (31).

The scale φ(t) and wavelet ψ(t) functions associated to the FIR filters {g n } and {h n } have compact support, thus offering the time resolution functionality. The simplest scale function that satisfies (25) is the characteristic function of the interval I = [0, 1), that corresponds to the Haar's scale function:

In this case (Haar MRA), the associated Haar scale filter is given by We can demonstrate that [17] :

and that

According to (42) and (43), we can obtain the coefficients u j (n) and w j (n) from the scale coefficients u j−1 (m) by means of decimation operation of the sequence {u j−1 (m)} by a factor of 2. The decimation consists in cascading a low-pass filter g(−m) (with a transfer functionḠ(z) = G(1/z) and frequency response G * (f )) or a high-pass h(−m) (with transfer functionH(z) = H(1/z) and frequency response H * (f )) with a compressor (or decimator) by a factor of 2. Decimate a signal by a factor D is the same as to reduce its sampling rate by D times.

The MRA is implemented by a low-pass and high-pass analysis filter banks G * (f ) and H * (f ) adequately positioned for separating the scale and wavelet coefficients sequences. This is known in the literature as the pyramid algorithm presented by Mallat [17] . Later, it is possible to rebuild the original signal using dual QMF reconstruction filter banks, low-pass G(f ) and high-pass H(f ).

It is important to emphasize that the pyramid algorithm's complexity is O(N ) (assuming we want to evaluate the DWT of N samples), while the direct evaluation of the DWT (that involves matrices multiplication) is O(N 2 ) [4] . Figure 8 shows the QMF analysis filter banks G * (f ) (low-pass) and H * (f ) (high-pass) with decimation (downsampling) by a factor of 2. Figure 9 shows the QMF reconstruction filter banks with interpolation (upsampling) by a factor of 2.

Observe that are used dual low-pass and high-pass filters, G(f ) and H(f ). Figure 10 presents the flow diagram that shows the initial projection of a signal x(t) on V 0 followed by the decomposition in W 1 , W 2 and V 2 . Figure 10 : Flow diagram that shows the initial projection of a signal x(t) on V 0 followed by the decomposition in W 1 , W 2 and V 2 . Figure 11 shows the flow diagram that illustrates the approximate synthesis of x(t) from W 1 , W 2 and V 2 . Figure 12 presents a block diagram that shows that the DWT works as a sub-bands codification scheme. The spectrum U 0 (f ) of the signal u 0 (n) is subdivided in three frequency bands (that cover two octaves): 0 ≤ f < 1/8, 1/8 ≤ f < 1/4 and 1/4 ≤ f ≤ 1/2. Figure 11 : Flow diagram that illustrates the approximate synthesis of x(t) from W 1 , W 2 and V 2 . Figure 12 : The spectrum U 0 (f ) of the signal u 0 (n) is subdivided in three frequency bands (that cover two octaves):

In section 2, we described how the DWT can be applied to a signal x(t). However, the purpose of this study requires that we think of x(t) as a realization of a stochastic process {x(t)}, so that we can be able to realize a statistical assessment of the change point. The reader interested in the definition of a stochastic process can consult the reference [24] for more information. Figure 13 shows the historical series in Brazil from January 2015 to July 2020 (67 samples). Figure 13 : Historical series in Brazil from January 2015 to July 2020.

In the paper [24] , we showed that the signal of Fig. 13 obeys a model of the type [25] [p. 58]

where x(t) are the observations, y(t) is a stationary process, and µ(t) denotes a linear trend given by the regression model

in which β 0 and β 1 are the intercept and the slope parameters. We also estimated an autoregressive model of order 11 (AR(11)) for y(t). An AR (11) is said to exhibit short-range dependence, as the Power Spectral Density (PSD) of the signal does not have 1/f behavior for frequencies near to zero (Long Range Dependence (LRD)).

In [24] , we conclude that there are no change points in the slope of the historical series and in the mean of the series that corresponds to the first difference of the historical series. Table 1 shows the estimated coeffcientes for (45) and its p-values. Figure 14 shows the historical series with the superimposed linear regression model. As in [24] , the first step in wavelet analysis requires that the deterministic trend of the historical series be removed, in order to perform the time-frequency domain analysis [25] [p. 58]. There are two alternatives: remove the line estimated by the regression or take the first difference in the series. As our goal is to coerce the data to (a possible) stationarity, then differencing may be more appropriate [25] [p. 61]. The first two samples of the historical series were discarded so that the series corresponding to the first difference has 64 points, that is, 2 6 points, which is appropriate for a wavelet analysis with the DWT.

The first difference can be denote as Figure 15 shows the series r(t) (we also demeaned it). Now, as far as a qualitative analysis is concerned, we have to focus our wavelet (graphical) analysis on the most refined scales of the time-frequency doman, where possible changes in variance can be localized. This is the intuition behind 

where we assume that r(t) has N j nonboundary wavelet coefficients w j (L j ), w j (L j + 1), . . . , w j (N j − 1), in which Here is a noteworthy example of an important advantage of wavelet analysis over Fourier. Should the same signal had been analyzed by the FT, we would not have been able to detect the instant of the change point, whereas it is clearly observable here. Now, let us compare the plot of the raw data along with the level-one wavelet detail signal (D1) in Fig. 18 (LA(8) MRA), which indicates the existence of a change in variance around t = 22 (2017/jan), see the red arrow. Remember that the detail signal D1 is located in the frequency band 1 4 , 1 2 . Fig. 19 illustrates the Haar MRA. Note that we can also see change in variance around t = 22 (2017/jan) in the level-one wavelet detail signal. Now we can confirme our qualitative analysis using the function wvarchg() (MATLAB© R2015a Wavelet Toolbox), which calculates the optimal positioning and (potentially) number of changepoints. As the signal r(t) is SRD, there is no need to be concerned with the question of sampling of the variance of the wavelet coefficients, which is problematic for long memory process [4, p. 380 ].

The test confirms that there is a change point in variance in t = 22.

In this paper, we presented an overview of the Continuos and Discrete Wavelet Transforms, and a wavelet analysis of the historical series of the total number of deaths per month in Brazil since 2015. Our preliminary results indicate that:

• the variance of the signal changed in january of 2017; however, this change point occurred before the outbreak of the COVID-19 pandemic; and • there is no evidence that COVID-19 provoked a change in the stochastic process that generates the historical series, as there are no change points in this signal after the outbreak of COVID-19 in Brazil, which happened in the beginning of march, 2020.

There is no doubt that COVID-19 has caused the deaths of many people around the world. Many who survived will have to live with sequelae in the brain, kidneys, lungs and heart. In addition, the world population is being subjected to a high degree of stress, due to fear of the new coronavirus, economic/financial problems, etc.

However, our results suggest that COVID-19 did not cause any change point in the brazilian series of the total number of deaths per month in Brazil since 2015 so far. What would be a plausible explanation for this strange result? Answering this question requires a multidisciplinary approach, as it involves several areas of knowledge such as medicine, signal processing, statistics etc.

Nevertheless, researchers in the area must continue to monitor the behavior of the historical series in Brazil.

In future work, we suggest the application of the Maximal Overlap Discrete Wavelet Transform (MODWT) and the Discrete Wavelet Packet Transform (DWPT) [4] to the brazilian historical series. We would also like to analyze the series from other countries affected by COVID-19.

Brazilian Federal Ministry of Health (Ministério da Saúde do Brasil). COVID-19 no Brasil

World Health Organization. Coronavirus disease (COVID-19) pandemic

Wavelet Methods for Time Series Analysis

Parametric statistical change point analysis

A wavelet-based approach for detecting changes in second order structure within nonstationary time series

Civil Registry Offices of Brazil (Cartórios de Registro Civil do Brasil). Transparency Portal (Portal da Transparência

The R Project for Statistical Computing

Code and Data

Theory of communication

A Friendly Guide to Wavelets

An Introduction to Wavelets and Other Filtering Methods in Finance and Economics

Orthonormal bases of compactly supported wavelets

Decomposition of hardy functions into square integrable wavelets of constant shape

A theory for multiresolution signal decomposition: The wavelet representation

Multiresolution approximations and wavelet orthonormal bases of l 2 (R)

Multifrequency channel decompositions of images and wavelet models

Meaningful MRA initialization for discrete time series

A Wavelet Tour of Signal Processing

Wavelet analysis of long-range dependent traffic

An exploratory time series analysis of total deaths per month in Brazil since

Time Series Analysis and Its Applications with R Examples