title: An exploratory time series analysis of total deaths per month in Brazil since 2015
authors: Alexandre Barbosa de Lima
text: 
  Tabular Structure Detection from Document Images for Resource Constrained Devices Using A Row Based Similarity Measure
      Ziheng Wang
      December 30, 2023
  ======================================================================================================================
  
  
  
  
  In this article, we investigate the historical series of the total number of deaths per month in Brazil since 2015 using time series analysis techniques, in order to assess whether the COVID-19 pandemic caused any change in the series' generating mechanism. The results obtained so far indicate that there was no statistical significant impact.
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Brazil is a country of continental dimensions, with a territorial extension of 8,510,295.914 Km^2, and an estimated population of 211,959,316 people, which is distributed in a federation constituted by 26 States and the Federal District (the country has 5,570 cities) <cit.>.
  Until the time of this writing, the Brazilian Federal Ministry of Health has recorded 3,622,861 Sars-CoV-2 case reports and 115,309 deaths <cit.> caused by the COVID-19 pandemic <cit.>.
  
  According to the Coronavirus Resource Center of the Johns Hopkins University (JHU) <cit.>, Brazil is the second country most affected by COVID-19 in the world, both in number of cases and in number of deaths.
  
  This article aims to investigate the historical series of the total number of deaths using time series analysis techniques. From now on, such series will be called historical series. More specifically, we want to assess the impact of COVID-19 on the evolution of the historical series. 
  
  The analysis was performed using the  software, version 4.0.2 <cit.>. The developed  code, as well as the database in Excel spreadsheet format, are available for public consultation  and auditing on GitHub <cit.>.
  
  The database used is that made available online by the Transparency Portal of the Civil Registry Offices of Brazil <cit.>, which consolidates the amount of birth, marriage and death certificates available in Brazil. Online data has been available since January 2015.
  Brazilian registries are regulated by the National Council of Justice (CNJ), which is a public institution headquartered in Brasília, Federal District, that aims to improve the work of the Brazilian judicial system, especially with regard to administrative and procedural control and transparency <cit.>. The president of the Brazilian Supreme Court also presides the CNJ.
  
  The remainder of the paper is organized as follows.  We review basic concepts of time series analysis  in Section <ref>. Section <ref> presents the technique used for modeling the historical series. This section also includes an exploratory data analysis and a  spectral analysis. Section <ref> summarizes the conclusions and highlights some topics for further investigation.
  
  
  
  
  
  § BACKGROUND
  
  
  
  
   §.§ Time Series Basic Concepts
  
  
  In broad terms, a time series consists in a set of numbers corresponding to the observation of a certain phenomenon. Figure <ref> shows the Nile River mimima for the years 622 AD to 1284 AD[Yearly minimal water levels of the Nile river measured at the Roda gauge near Cairo.] <cit.>. By nature, such numbers are realizations of random variables. In general, a collection of random variables, {_t}, indexed by t, is referred to as a stochastic process <cit.>. In this paper, t will be discrete and vary over the integers t=0,± 1, ± 2, …. 
  
  
  
  Box and Jenkis <cit.> introduced the class of stationary ARMA(p,q) (autoregressive moving average) models[In this work, we use the simplified notation _t to denote a discrete-time stochastic process {_t}.]
  
  
      _t - μ = ∑_j=1^pϕ_j (_t-j-μ) + _t - ∑_i=1^qθ_i _t-j,
  	
  
  where μ=E{_t} is the mean of _t, {ϕ_1, ϕ_2,…,ϕ_p} and {θ_1, θ_2,…,θ_q} are parameters of the model, and _t is a wide-sense stationary white noise process with zero mean and power σ^2, i. e., _t ∼ (0,σ^2). In a more compact form, we have
  
  
      ϕ(B) '_t = θ(B)_t,
  	
  
  where '_t = _t - μ, B is the backward shift operator (B_t = _t-1), ϕ(B) is the autoregressive (AR) operator of order p
  
  
      ϕ(B) = 1 - ϕ_1 B - ϕ_2 B^2 - … - ϕ_p B^p
  
  
  and θ(B) denotes the moving average (MA) operator of order q
  
  
      θ(B) = 1 - θ_1 B - θ_2 B^2 - … - θ_q B^q.
  	
  
  In the rest of this paper, we will assume μ=0 without loss of generality.
  
  The process _t can be viewed as the output of a digital filter (ARMA filter) whose input is _t, with system function 
  
  
      H(z) = θ(B)/ϕ(B),
  
  
  where H(z) denotes the z-transform of the impulse response h_t of the ARMA filter. An ARMA(p,q) process _t is said to be wide-sense stationary (or non-explosive)  if the poles of H(z) in (<ref>) lie inside the complex unit circle (|z|=1), and it is invertible if the zeros of H(z) in (<ref>) lie inside the unit circle. 
  The autocorrelation function (ACF) [We assume that the autocorrelation function is given by  ρ_h=γ_h/γ_0, where γ_h corresponds to the autocovariance of _t at lag h.] ρ_h  of an ARMA(p,q) process shows exponentially decay to zero, i. e., at lag h converges rapidly to zero as h→∞ (short memory property) <cit.>. 
  
  A random process _t is wide sense (or weakly) stationary if its mean is constant <cit.>
  
  
      E[_t] = μ_,
  
  
  and if its ACF depends only on the lag τ=t_2-t_1:
  
  
      ρ_(t_1,t_2) = ρ_(t_1,t_1 + τ) = ρ_(τ).
  
  
  In the literature, it is common to use the terms time series and stochastic process interchangeably <cit.>, <cit.>, <cit.>, <cit.>. From now on, we will only use the term time series. The context will indicate to the reader whether it is a process or a realization of a process.
  
  
  
   §.§ Time Series Modeling
  
  
  The modeling of a time series _t consists on estimating an invertible function h(.), called model of _t, such that
  
  
      _t = h(…,_t-2,_t-1,_t,_t+1,_t+2,…),
  
  
  in which _t∼Independent and Identically Distributed – IID and
  
  
      g(…,_t-2,_t-1,_t,_t+1,_t+2,…) = _t,
  
  
  in which g(.)=h^-1(.). The process _t is the innovation at instant t and represents the new information about the series that is obtained at instant t. 
  
  In practice, the adjusted model is causal, i. e., 
  
  
      _t = h(_t,_t-1,_t-2,…).
  
  
  The model construction methodology is based on the iterative cycle illustrated by Fig.  <ref> <cit.>:
  
  
  	
    (a) a general class of models is considered for the analysis (specification);
  	
    (b) there is the identification of a model, based on statistical criteria; 
  	
    (c) it follows the estimation phase, in which the model's parameters are obtained. In practice, it is important that the model is parsimonious[We say that a model is parsimonious when it uses few parameters. The use of an excessive number of parameters is undesirable because the uncertainty degree of the statistical inference procedure increases with the number of parameters.];  and
  	
    (d) at last, there is the diagnostic of the adjusted model by means of a statistical analysis of the residual series w_t (is _t compatible with a white noise process?)
  
  
  
  
  The process _t of (<ref>) is linear when it corresponds to the convolution of a process _t ∼IID and a deterministic sequence h_t <cit.>[p. 377] 
  
  
      _t    = h_t ⋆_t = ∞k=0∑h_k _t-k
         = _t + h_1_t-1 + h_2_t-2 + …
         = (1 + h_1B + h_2B^2 + …)_t 
         = H(B)_t
  
  
  in which the symbol ⋆ denotes the convolution operation and h_0=1. 
  
  Eq. (<ref>) is also known as the infinite order moving average (MA(∞)) representation <cit.>.
  
  As the models in practice are invertible, the model of _t  can be rewritten in an infinite order autoregressive  (AR(∞)) form:
  
  
      _t  = ∑_k=1^∞ g_k _t-k + _t.
  	
  
  An order p AR model satisfies the equation
  
  
      ϕ(B) _t = _t.
  
  
  in which ϕ(B) is an order p polynomial.
  
  
  
   §.§ Identifying AR Models
  
  
  In practice, the order p of an AR series is unknown and must be empirically specified. In this paper, we use an information criterion function <cit.> as will be explained below.
  
  The basic idea of an ARMA model selection criterion is to choose the orders k and l that minimize the quantity 
  
  
      P(k,l) = lnσ̂^2_k,l + (k+l)C(N)/N,
  
  
  in which σ̂^2_k,l is a residual variance estimate obtained by adjusting an ARMA(k,l) model to the N series observations, and C(N) is a function of the series size. 
  
  The quantity (k+l)C(N)/N is called penalty term and it increases when the number of parameters increases, while σ̂^2_k,l decreases.
  
  Akaike proposed the  information criterium <cit.>, <cit.>
  
  
      AIC(k,l) = lnσ̂^2_k,l + 2(k+l)/N,
  
  
  known as AIC, in which σ̂^2_k,l is the maximum likelihood estimator of  σ^2_ for an ARMA(k,l) model. 
  
  Upper bounds K and L for k and l must be specified. Eq. (<ref>) has to be evaluated for all possible (k,l) combinations with 0≤ k ≤ K 0≤ l ≤ L. In general, K and L are functions of N, for example, K=L=ln N.	
  
  For the case of AR(p) models, (<ref>) reduces to
  
  
      AIC(k) = lnσ̂^2_k + 2k/N,    k ≤ K.
  
  
  
  
   §.§ Estimation of AR Models
  
  
  Having identified the AR model's order p, we can go to the parameters estimation phase. The methods of moments, Least Squares and Maximum Likelihood may be used <cit.>, <cit.>. As, in general, the moments estimators are not good <cit.>, statistical packages as  and  use some Least Squares or Maximum Likelihood estimator.
  
  
  
   §.§ ARIMA model
  
  
  If a process which corresponds to the difference of order d=1,2,… of _t
  
  
      _t= (1-B)^d _t = Δ^d _t
  
  
  is stationary, then _t can be represented by an ARMA(p,q) model
  
  
      ϕ(B) _t = θ(B)_t.
  
  
  In this case, 
  
  
      ϕ(B)Δ^d _t = θ(B)_t
  
  
  is an ARIMA(p,d,q) model and we say that _t is  an “integral” of _t <cit.> because 
  
  
      _t= S^d _t.
  
  
  The ARIMA(p,d,q) model
  
  
      H(z)= θ(z)/ϕ(z)(1-z^-1)^d
  
  
  is marginally stable <cit.>, as it has d roots on the unit circle. Also, _t of (<ref>) is a homogeneous non-stationary process (meaning non-explosive) or having  unit roots <cit.>, <cit.>, <cit.>.
  
  Observe that <cit.>[p.139]:
  
  
  	
    (a) d=1 corresponds to homogeneous non-stationary series with respect to the level (they oscillate around a mean level during a certain time and then jump to another temporary level);
  	
    (b) d=2 corresponds to homogeneous non-stationary series with respect to the trend (they oscillate along a direction for a certain time and then change to another temporary direction).
  
  
  The ARIMA model (<ref>) may be represented in three ways:
  
  
  	
    (a) ARMA(p+d,q) (similar to Eq. (<ref>)) 
  	
      (t) = p+dk=1∑φ_k (t-k) + (t) - qk=1∑θ_k (t-k);
  
  	
    (b) AR(∞) (inverted format), given by (<ref>) or
  	
    (c) MA(∞), according to (<ref>).
  
  
  
  
   §.§ Random walk
  
  
  Consider the model _t ∼ I(1) 
  
  
      _t = _t-1 + _t,
  
  
  in which _t is a stationary process. If we assume the initial condition y_0,  (<ref>) can be rewritten as an integrated sum
  
  
      _t = _0 + ∑_j=1^t_j.
  
  
  The integrated sum ∑_j=1^t_j is called stochastic trend and it is denoted by  TS_t. Observe that
  
  
      TS_t = TS_t-1 + _t,
  
  
  in which TS_0=0.
  
  If _t ∼𝒩(0,σ^2_) in (<ref>), then _t is known as random walk. 
  
  Including a constant in the right side of  (<ref>), we have a random walk with drift,
  
  
      _t = θ_0 + _t-1 + _t.
  
  
  Given the initial condition y_0, we can write
  
  
      _t    = y_0 + θ_0 t + ∑_j=1^t_j
         = TD_t + TS_t
  
  
  The mean, variance, autocovariance and ACF of _t are given by <cit.>
  
  
      μ_t   = y_0 + t θ_0 
      σ^2(t)               = t σ^2_
      
      C_k(t)               = (t-k)σ^2_
      ρ_k(t)            = t-k/t.
  
  
  Observe that ρ_k(t)≈ 1 when t>>k and the literature states that the random walk has “strong memory” <cit.>.
  
  The random walk's Sample Autocorrelation Function (SACF) decays linearly for large  lags. 
  
  
  
   §.§ Spectral Analysis of Random Signals
  
  
  Spectral analysis is a well-established research area <cit.>. However, the estimation of the power spectrum of a signal is not a trivial matter. There are two classes of spectral analysis techniques currently in use:  parametric (or model-based) and nonparametric analysis. Both methods are used in this work.
  
  The fundamental idea of parametric spectral analysis is fairly simple. The parametric approach assumes that the signal satisfies a generating model, such as an AR(p), with known functional form and then proceed by estimating the parameters in the assumed model <cit.>, <cit.>. The most widely used form of parametric Power Spectral Density (PSD) estimation uses an AR(p) model <cit.>.
  
  Let us now consider the nonparametric (or classical) method. Then, the estimation of PSD of a time series _t can be made using periodogram methods based on the Discrete Fourier Transform (DFT), which can be efficiently calculated by a Fast Fourier Transform (FFT) Algorithm. 
  
  In the sequence we present the periodogram method.
  
  Consider a time series _t with N values or samples, i. e., _t=0 outside the time interval 0 ≤ t ≤ N-1. In some cases of interest, we consider that _t has a size N, even if its actual size is M ≤ N (in such cases the series _t must be completed with (N-M) zeros (zero padding).
  
  The equation (<ref>) is the DFT of _t: 
  
  
      X[k]= {[ N-1t=0∑_t e^-j 2πk/N t,               0≤ k≤ N-1;                      0,               otherwise ].
  
  
  The PSD of _t can be estimated by calculating the periodogram, given by
  
  
      P_x(f_k) = |X[k]|^2/N
  
  
  where X[k] denotes the DFT of _t, and 
  
  
      f_k=0,(1/N), …, (k/N), …, (N-1/N)
  
  
  The periodogram is an asymptotically unbiased spectral estimator of the PSD of a random signal <cit.>, <cit.>. Its main problem lies in its large variance. In other words, the periodogram is inconsistent (i. e., the dispersion of the estimates is independent of N). This motivates the use of a “refined periodogram method”, like the Daniell method.
  
  It is possible to show that the periodogram values P_x(f_k) are asymptotically uncorrelated random variables <cit.>. Thus we may reduce its large variance by weight averaging the periodogram over small intervals centered on the current frequency f_k. The practical form of the Daniell estimate can be performed using the FFT. This work uses a Daniell kernel for nonparametric PSD estimation. For further details, please refer to <cit.> and <cit.>.
  
  
  
  § EXPERIMENTAL RESULTS
  
  
  
  
   §.§ Time Series Regression Modeling
  
  
  Figure <ref> shows the historical series in Brazil from January 2015 to July 2020 (67 samples).
  
  
  
  The series in the Fig. <ref> strongly suggests that we can specify a model of the type <cit.>[p. 58]
  
  
      _t  = μ_t + _t
  
  
  where _t are the observations, _t is a stationary process, and μ_t denotes a linear trend given by the regression model
  
  
      μ_t = β_0 + β_1 t
  
  
  in which β_0 and β_1 are the intercept and the slope parameters.
  
  Table <ref> shows the estimated coeffcientes for (<ref>) and its p-values, which are negligible. 
  The goodness of fit is summarized by the R^2 of the regression <cit.>[p. 169]. Note that the  R^2 statistics given by Table <ref> indicate that the proposed model for μ_t explain approximately 75.8% of the total variability of the data. The great value of the F-statistic and the neglibible model p-value show that we can not reject the null hypothesis of linear regression. Thus, we can consider a linear model for (<ref>) to be statistically significant given the statistical significance level of 0.01.  Figure <ref> shows the time series with the superimposed linear regression model.
  
  
  
  
  
  
  
  
  To verify that the model (<ref>) is appropriate, it is also necessary to investigate the residuals. This is what is called residual analysis.
  
  The residuals of (<ref> ) correspond to the discrepancies between the observed values (μ) and the values adjusted (μ̂) by the model. The i-th residual is given by
  
  
      ê_i = μ_i - μ̂_i.
  
  
  Figure <ref> shows  the residuals vs fitted model and the Quantile-Quantile plot (QQ-plot) of the residuals, which suggests that they are normally distributed.
  
  
  
  Table <ref> shows a residual diagnostics using the Jarque-Bera and Shapiro-Wilks tests <cit.>[p.61] for the null hypothesis of normality of the residuals. The null of normality is not rejected using either tests.
  
  
  
  At this point, we can conclude that:
  
  
  	
    * there is no statistical evidence that COVID-19 affected the deterministic linear trend of the historical series, i. e., on average, the monthly growth in the number of deaths, which is approximately 675 deaths/month, did not change since the first recorded death in Brazil on March 16, 2020; and
  	
    * there is no change point in the deterministic linear trend of the historical series.
  
  
  
  
   §.§ Exploratory Data Analysys
  
  
  The first step in exploratory analysis is to remove the deterministic trend of Eq.(<ref>) <cit.>[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 <cit.>[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 facilitates the spectral analysis with FFT.
  
  The first difference can be denote as
  
  
      Δ_t = _t - _t-1 = _t.
  
  
  Figure <ref> shows the series _t (we also demeaned the series) and its SACF. 
  
  Figure <ref> shows the histogram of _t with a superimposed normal distribution and its Q-Q plot, which suggests that _t  follows a normal distribution.
  
  Figure <ref> shows the smoothed periodogram using the Daniell method and the PSD for the estimated model of _t, which is an AR(11).
  
  We used the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test <cit.> for the null hypothesis that _t is level stationary, i. e., that the series is I(0). We obtained a p-value of 0.1, which means that one can not reject that _t is level stationary.
  
  Finally, but not least, an informal analysis of Fig. <ref> suggests that the variance of _t has a change point around t=22 (2017/jan). However, if this really happened, which we cannot guarantee with the techniques employed in this article, ocurred long before the outbreak of COVID-19 in Brazil. This possible change point motivates a time-frequency domain analysis using wavelets, as the localized nature of wavelet coefficients allows one to analyze the evolution of the series variance over time <cit.>. This will be investigated in future work.
  
  
  
  
  
  
  
  
  
  
  § CONCLUSIONS AND FUTURE WORK
  
  In this paper, we presented an exploratory time series analysis of the historical series of the total number of deaths per month in Brazil since 2015. Our preliminary results indicate that: 
  
  
  	
    * there is no statistical evidence that COVID-19 affected the deterministic linear trend of the historical series, i. e., on average, the monthly growth in the number of deaths, which is approximately 675 deaths/month, did not change since the first recorded death in Brazil on March 16, 2020; 
  	
    * there is no change point in the deterministic linear trend of the historical series; 
  	
    * there is significant statistical evidence that the first difference time series is stationary; and
  	
    * there is no statistical evidence that COVID-19 provoked a change in the stochastic process that generates the time series under analysis[The random process is the series' generating mechanism.].	
  
  
  These results are thought provoking and not intuitive. COVID-19 has caused many deaths around the world. This is an indisputable fact. However, our results suggest that this disease does not have so far an additive effect on the total number of deaths per month in Brazil. What would be a plausible explanation for this strange result?
  
  In any case, further research should be carried out to confirm the results obtained.
  
  In future work, we will analyze the historical series using wavelet methods. 
  
  
  
  IEEEtran
  
  
  