There is a growing interest in the last years to deal with data that is only partially registered or underreported in the time series literature. This phenomenon is very common in many fields, and has been previously explored by different approaches in epidemiology, social and biomedical research among many other contexts [1] [2] [3] [4] [5] . The sources and underlying mechanisms that cause the underreporting might differ depending on the particular data. Some authors consider a situation where the registry is updated with time and therefore the underreporting issue is mitigated 6 . That leads to temporary underreporting while this work is focused on permanent underreporting, where the registered data are never updated in order to become more accurate. From the methodological point of view, several alternatives have been explored, from Markov chain Monte-Carlo based methods 5 to recent discrete time series approaches 7, 8 and several attempts to estimate the degree of underreporting in several contexts have been done 9 , although there is a lack of models incorporating continuous time series structures and handling underreporting.
One of the fields where the interest in addressing the underreporting issues is higher is the epidemiology of infectious diseases. In the last few years, many approaches to deal with underreported data have been suggested with a growing level of sophistication from the usage of multiplication factors 10 to several Markov-based models 11, 12 or even spatio-temporal modelling 13 . Even a new R 14 package able to fitting endemic-epidemic models based on approximative maximum likelihood to underreported count data has been recently published 15 .
Human papillomavirus (HPV) is one of the most prevalent sexually transmitted infections. It is so common that nearly all sexually active people have it at some point in their lives, according to the information provided by the Centre of Control of Diseases (CDC) 16 . Generally, the infection disappears on its own without inducing any health problem, but in some cases it can produce an abnormal growth of cells on the surface of the cervix that could potentially lead to cervical cancer. HPV infection is also related to several cancers (cervical, vulva, vagina, penis, anus, . . . ) and other diseases like genital warts (GW). The fact that most cases of HPV infection are asymptomatic causes that public health registries might be potentially underestimating its incidence. The underreporting phenomenon in HPV data from the discrete time series point of view has been recently studied 7 .
There is an enormous global concern around 2019novel coronavirus (COVID-19) infection in the last few months, leading the World Health Organization (WHO) to declare public health emergency 17 . As the symptoms of this infection can be easily confused with those of similar diseases like Middle East Respiratory Syndrome Coronavirus (MERS-CoV) or Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV), its incidence has been notably underreported, especially at the beginning of the outbreak in Wuhan (Hubei province, China) by December 2019.
The characteristics of the proposed model and the methodology to estimate the parameters are detailed in Section Model definition. The application of the model on the HPV incidence data in Girona (Catalunya, Spain) in the period 2010-2014 is illustrated in Section Example: HPV infection incidence and on the COVID-19 incidence in Heilongjiang (China) is shown in Section Example: Coronavirus COVID-19 incidence in the region of Heilongjiang. Finally the results of a thorough simulation study are summarised in Section Simulation study.
Consider an unobservable process with an AutoRegressive Moving Average (ARM A(p, r)) structure defined by
where ǫ t is a Gaussian white noise process with ǫ t ∼ N (µ ǫ , σ 2 ǫ ). In our setting, this process cannot be directly observed, and all we can see is a part of it, expressed as
The interpretation of the parameters in Eq. (2) is straightforward: q is the overall intensity of misreporting (if 0 < q < 1 the observed process Y t would be underreported while if q > 1 the observed process Y t would be overreported). The parameter ω can be interpreted as the overall frequency of misreporting (proportion of misreported observations).
Consider that the unobserved process X t follows an ARM A(p, r) model as defined in Eq. (1). As can be seen in Appendix 1, the observed process has mean
The autocorrelation function of the observed process can be written in terms of the features of the hidden process X t as
where ρ X is the autocorrelation function of the unobserved process X t .
A situation of particular interest is the case ω = 1, meaning that all the observations might be underreported and that a simpler model for Y t excluding the parameter ω might be suitable
In this case, however, the observed process Y t would be a non-identifiable ARM A(p, r) as parameter q could not be estimated on the basis of the methodology described in the following section.
The likelihood function of the observed process Y t is not easily computable but the parameters of the model can be estimated by means of an iterative algorithm based on its marginal distribution, using R packages mixtools 18 and forecast 19, 20 . The main steps are described in detail below:
(i) Following Eq. (2), the observed process Y t can be written as
where Z t is an indicator of the underreported observations, following a Bernoulli distribution with probability of success ω (Z t ∼ Bern(ω)), its marginal distribution is a mixture of two normal random variables N µ, σ 2 and N q · µ, q 2 · σ 2 respectively, where
. This fact can be used to obtain initial estimates for q and ω. Using the E-M algorithm (specifically on the Estep), the package mixtools calculates the posterior probabilities (conditional on the data and the obtained estimates) of each observation to come from one of these two normals.
(ii) Using the indicatorẐ t obtained in the previous step, the series is divided in two: One including the underreported observations (treating the nonunderreported values as missing data) and another with the non underreported observations (treating the underreported values as missing data). An ARIM A model is fitted to each of these two series and a newq is obtained by dividing the fitted means. (iii) A mixture of two normals is fitted to the observed series Y t with mean and standard deviation fixed to the corresponding values obtained from the previous step, and a new ω is estimated. (iv) Steps (ii) and (iii) are repeated until the quadratic distance between two consecutive iterations
Once the parameter estimates are stable according to the previous criterion, the underlying process X t is reconstructed asX
and an ARIM A model is fitted to the reconstructed process to obtainα j , j = 1, . . . , p,θ k , k = 1, . . . , r andσ ǫ 2 .
To account for potential trends or seasonal behaviour, covariates can be included in the described estimation process. Additionally, a parametric bootstrap procedure with 1000 replicates is used to estimate standard errors and build confidence intervals based on the percentiles of the distribution of the estimates. All the R codes used to estimate the parameters and build the confidence intervals are available as supplementary material.
The results of the proposed methodology over a comprehensive simulation study and an application on two real data sets are shown in this Section.
A thorough simulation study has been conducted to ensure that the model behaves as expected in several situations, including AR(1), M A(1) and ARM A(1, 1) structures for the hidden process X t with values for the parameters α, θ, q and ω ranging from 0.1 to 0.9 for each parameter. For each autocorrelation structure and parameters combination, a random sample of size n = 1000 has been generated using the function arima.sim from R package forecast 19, 20 . Different sample sizes (n = 50, 100, 500) have also been analysed to study the impact of sample size on accuracy and the results are reported in the Supplementary Material. The complete tables with the estimates corresponding to each combination of parameters can be found in Table S4 to S6 n = 50 and AR(1), M A(1) and ARM A(1, 1) respectively, Tables S7 to S9 for n = 100 and Tables S10 to S12 for n = 500, and the performance is summarised in Tables S13 to S15 for n = 50, 100 and 500 respectively. Absolute average bias is similar regardless of the sample size, while average interval lengths (AIL) are higher and interval coverages are poorer (around 75% for n = 50) for lower sample sizes as could be expected. Several bootstrap sizes (b = 20, 50, 100) were also considered and the difference between them were negligible, so only results corresponding to b = 100 bootstrap replicates are reported. The average absolute bias, interval coverage and 95% confidence interval length are reported in Table 1 . These values are averaged over all combinations of parameters. Additionally, standard AR(1), M A(1) and ARM A(1, 1) models were fitted to the same simulated series without accounting for their underreporting structure. The complete table with the estimates corresponding to each combination of parameters is available as Supplementary Material (Tables S1, S2 and S3 for AR(1), M A(1) and ARM A(1, 1) structures respectively).
It is clear from Table 1 that ignoring the underreported nature of data (labeled as Standard models in the table) leads to highly biased estimates with extremely low coverage rates, even with larger average interval lengths. This is especially relevant when the intensity or frequency of underreported observations is high.
The series of weekly cases of HPV infection in Girona in the period 2010-2014 was previously analysed as a discrete INAR(1) hidden Markov process 7 . In a similar way, we aim to analyse the corresponding series of incidence, and an AR process of order 1 seems to be adequate (see Figure 1 ). According to Eq. (3), the autocorrelation function of the observed process Y t when the hidden process X t has an AR(1) structure takes the form ρ Y (k) = c · α k , where c = c(α, µ ǫ , σ 2 ǫ , ω, q) =
In particular, in this case we can write log(ρ Y (k)) = log(c) + k · log(α), so a statistically significant intercept of this linear regression model (estimating the parameters by ordinary least squares method) could be interpreted as an evidence of underreporting, as in this case (p − value = 0.0014). It is clear from Figure 1 that the estimated regression line does not cross the origin, so the behavior of the observed process is consistent with an underlying underreported AR(1) process. By means of the estimation method described in Section Estimation, it can be seen that the estimated model for the hidden process is X t = 0.104 · X t−1 + ǫ t , being the observed process Y t , Y t = X t with probability 0.202 0.236 · X t with probability 0.798 (5) The estimated parameters are reported in Table 2 . These results are highly consistent with those previously reported in the literature for the number of HPV cases obtained through a discrete time series approach 7 and can be interpreted in a very straightforward way. Moreover, this new methodology can be used to model the incidence of the disease instead of the number of cases, accounting for potential changes in the underlying population.
The estimated intensity of underreporting isq = 0.236, with 95% confidence interval (0.150, 0.322). The registered and estimated evolution of HPV incidence within the study period (2010-2014) can be seen in Figure 2 . 
COVID-19 is a betacoronavirus that affects the lower respiratory tract and often manifests as pneumonia in humans which was identified as the causative agent of an unprecedented outbreak of pneumonia in Wuhan City, Hubei province in China starting in December 2019 17 . Considering that many cases run without developing symptoms beyond those of MERS-CoV, SARS-CoV or pneumonia due to other causes, it is reasonable to assume that the incidence of this disease has been underregistered, especially at the beginning of the outbreak 21 .
Heilongjiang is a province in north-east China. Although in general the behavior of this kind of diseases is far from being stationary, this province is far enough from the focus in Hubei province (south central China) so the incidence is much lower and less explosive, and in the study period of time (2020/01/22-2020/02/26) it can be considered stationary, as can be seen in Figure 4 , which shows registered and estimated evolution of COVID-19 incidence within the considered period of time.
A disease with a similar behavior as the MERS-CoV has been modeled as an ARM A(3, 1) 22 , so we checked this model and similar ones. Probably due to the shortness of the available data this autoregressive structure was not observed and in our case the best performing model was an M A(1) (AIC of -151.31 against -148.82 for the ARM A(3, 1)), consistently with the residuals profile shown in Figure 3 , obtained from fitting an M A(1) model to the most likely process X t reconstructed following step (v) in Section Estimation.
By means of the estimation method described in Section Estimation, it can be seen that the estimated model for the hidden process is X t = 0.481 · ǫ t−1 + ǫ t , being the observed process Y t , Y t = X t with probability 0.507 0.195 · X t with probability 0.493 (6) The estimated parameters are reported in Table 3 .
In biomedical and epidemiological research, the usage of disease registries in order to analyse the impact and incidence of health issues is very common. However, the accuracy and data quality of such registries is in many cases at least doubtful. This is the case, for instance, for rare diseases 23 or health issues that clear asymptomatically in most cases like HPV infection. In the case of HPV incidence in Girona in the period 2010-2014, the registered weekly average is 0.17 cases per 100,000 individuals, while the reconstructed process has a weekly average of 0.51 cases per 100,000 individuals, showing that only 33% of the real incidence is recorded by the public health system. This result is very consistent with that of 7 , claiming that only 38% of the HPV cases were registered in the same area and period of time.
The Heilongjiang region COVID-19 data reveal that in average about 66% of the actual incidence in the period 2020/01/22-2020/02/26 was reported. The unavailable data estimated by the proposed methodology are crucial to provide public health decision-makers with reliable information, which can also be used to improve the accuracy of dynamic models aimed to estimate the spread of the disease 21 .
The concerns around accuracy of registered data have recently led to the publication of recommendations to improve data collection to ensure accuracy of registries (see for instance 24, 25 ). Nonetheless, these recommendations are very recent and may be difficult for the public health services of many countries to fully implement them, due to operational or structural issues.
The proposed methodology is able to deal with underreported (or overreported) data in a very natural and straightforward way, estimating its intensity and frequency on a continuous time series, and allowing to reconstruct the most likely unobserved process. It is also flexible enough to handle covariates straightforwardly, and therefore it is simple to introduce trends or seasonality if necessary, so it can be useful in many contexts, where these issues might arise.
The simulation study shows that the proposed methodology behaves as expected and that the parameters used in the simulations, under different autocorrelation structures, are properly recovered, regardless of the intensity and frequency of the underreporting issues. and that X t and Z t are independent,
To compute the variance of the observed process Y t it is important to observe that E(
we only need to compute E((1 − Z t + qZ t ) 2 ):
and so E( 
The covariance of the hidden process X t can be written as 26
where ψ(z) = ∞ j=0 ψ j z j = θ(z) α(z) for | z |≤ 1, with θ(z) = 1 + θ 1 · z + . . . + θ r · z r and α(z) = 1 + α 1 · z + . . . + α p · z p . Considering Z t ∼ Bern(ω) we can write E(Y t Y t+k ) = E(X t · (1 − Z t )X t+k · (1 − Z t+k ))+ +E(X t · (1 − Z t )qX t+k · Z t+k )+ +E(X t+k · (1 − Z t+k )qX t · Z t )+ +E(qX t qX t+k · Z t Z t+k ) (11) As Z t and Z t+k are independent of {X j } and both have expectation ω,
(12) Therefore, E(Y t Y t+k ) = (1 − ω + q · ω) 2 E(X t X t+k ). From here, the covariance of the observed process Y t is
Finally, the autocorrelation of Y t can be expressed as ρ Y (k) = γY (k) V (Yt) = = γX (k)·(1−ω+q·ω) 2 (V (Xt)+E(Xt) 2 )·(1+ω·(q 2 −1))−E(Xt) 2 ·(1−ω+q·ω) 2 = = V (Xt)·ρX (k)·(1−ω+q·ω) 2 (V (Xt)+E(Xt) 2 )·(1+ω·(q 2 −1))−E(Xt) 2 ·(1−ω+q·ω) 2 ,
where ρ X (k) is the autocorrelation function of the hidden process X t .