key: cord-0809757-xeimo0go authors: Gourieroux, C.; Jasiak, J. title: Time varying Markov process with partially observed aggregate data: An application to coronavirus() date: 2020-11-28 journal: J Econom DOI: 10.1016/j.jeconom.2020.09.007 sha: 40760da5242de36bd61660c496d49c8475b89bdd doc_id: 809757 cord_uid: xeimo0go A major difficulty in the analysis of Covid-19 transmission is that many infected individuals are asymptomatic. For this reason, the total counts of infected individuals and of recovered immunized individuals are unknown, especially during the early phase of the epidemic. In this paper, we consider a parametric time varying Markov process of Coronavirus transmission and show how to estimate the model parameters and approximate the unobserved counts from daily data on infected and detected individuals and the total daily death counts. This model-based approach is illustrated in an application to French data, performed on April 6, 2020. The aim of this paper is to address the problem of partial observability, encountered recently in epidemiological research on Covid-19. More specifically, some individuals are infected and asymptomatic. Therefore, they remain undetected and unrecorded, especially during the early phase of the epidemic 1 . As a consequence, the total count of recovered and immunized individuals is unknown, as only the number of recovered detected individuals is available. This problem of partial observability of counts renders difficult the estimation of an epidemiological SIRD (Susceptible, Infected, Recovered, Deceased) model, extended to disentangle the infected and undetected from the infected and detected individuals. Moreover, such substantial undocumented infection can facilitate fast transmission of the virus (Li et al.(2020) ). The unknown total counts of infected individuals can be approximated by sampling the population daily and performing serological tests on the sampled individuals to estimate the rates of infected undetected and recovered individuals. However, it takes time to validate and produce reliable serological tests for Covid-19. Moreover, regularly performed sampling can be costly, especially in terms of time of health care providers. The alternative method, proposed in this paper, is purely model-based. Loosely speaking, under the standard extended SIRD model, the evolution of death rates might be different, depending on whether all infected individuals are detected or not. This implied difference will allow us for a model-based estimation of the proportions of infected undetected individuals (resp. recovered immunized) [see, Verity et al.(2019) for pure model based estimation of coronavirus infection, Manski, Molinari (2020) for set estimation of the infection rate]. This paper discusses the general case of time varying Markov processes when aggregate counts are partially observed. It is organized as follows. Section 2 describes the latent model of qualitative individual histories. These histories follow a time varying Markov process with transition probabilities that can depend on latent counts and unknown parameters. The observations are functions of the frequencies of individual states (called compartments in epidemiology), although not all of those frequencies are observed, in general. More specifically, only some states can be observed and/or a sum of frequencies Let f (t) denote the cross sectional frequency, i.e. the sample counterpart of p(t). It follows from the standard limit theorem that: Proposition 1: Under Assumption A1, the frequencies f (t) are consistent of p(t) and asymptotically normal for large N . Their variance-covariance matrix is given in Appendix 1. This specification of the transition matrix includes the homogeneous Markov chain, when there is no effect of lagged p(t−1). It also includes the standard contagion SIR-type models used in epidemiology [see, McKendrick (1926) , Kermack, McKendrick (1927) for early articles on SI and SIR models in the literature, Hethcote (2000), Brauer et al.(2001) , Vinnicky, White (2010) for general presentations of epidemiological models, Allen (1994) for their discrete time counterparts, Gourieroux, Jasiak (2020) for an overview, and also examples given below]. As vectors p(t) change over time, stationarity is not assumed. In practice, the individual histories, or the counts of flows between the states 2 may not be observed, while cross-sectional frequencies are generally available. These can be the frequencies f (t), t = 1, ..., T , or aggregates of such frequencies. Assumption A2: The observations are: t = Af (t), t = 1, .., T , where A is a K × J state aggregation matrix, that is a matrix with rows containing zeros and ones. The aggregation matrix is known and of full rank K. Although, in general, the process of aggregate counts: f 1 (t) + f 2 (t), f 3 (t) + f 4 (t) may not be Markov, it is important to consider the special case when it is, and then explore the possibility of identifying the parameters of the regional, i.e. disaggregated dynamics. This is the objective of the percolation theory [see, Garet, Marchand (2006) for a detailed analysis of competing contagion sources]. Under Assumption A1, we can use the Bayes' theorem to link the marginal theoretical probabilities p(t) to the transition probabilities as follows: for the sequence of cross-sectional distributions. These equations will be used as the estimating equations in the asymptotic least squares estimation method outlined below 4 . In our framework, the parameter of interest includes θ as well as the (equilibrium) sequence of vectors p(t). They can be jointly estimated from the following optimization: where ||.|| denotes an Euclidean norm. This estimation is constrained to account for the positivity and unit mass restrictions on the p(t)'s, and for potentially other restrictions on parameter θ (see, Section 5.3). The estimation method depends on the selected norm, such as ||p|| 2 = p p, or ||p|| 2 = p Ω −1 p, where Ω is a symmetric positive definite matrix, or a norm, which varies during the disease transmission depending on the precision of frequencies f (t) [see, Gourieroux, Jasiak (2020) ]. Proposition 2: If the constrained optimization given above has a unique solution, which is continuously differentiable with respect to Af (t) = t , t = 1, ..., T , then the estimator is asymptotically consistent, converges at rate 1/ √ N , and is asymptotically normally distributed. Proof: see Appendix 2. The expression of the asymptotic variance-covariance matrix is derived by a delta method from the asymptotic variance-covariance matrix of f (t) given in Appendix 1. If A = Id, that is, if all frequencies are observed, we obtain the case analysed in McRae (1977) . In the general framework, this optimization is not only used to estimate parameter θ, but also to approximate the unobserved marginal probabilities. For ease of exposition, let us consider A = Id. The constrained optimization (3.2) can be interpreted in a (pseudo) state space framework with the measurement equation: (3.4) and an assumption on the variance of u(t), depending on the selected Euclidean norm. This is a pseudo-state space representation, rather than an exact state space representation, as the errors u(t) are serially dependent [see, Appendix 1]. A Kalman filter 5 can be applied to the above pseudo-state space [for example, under the assumption of independent errors u(t) ∼ N (0, Σ)] to estimate numerically equation (3.2). However, the estimated elements of the variance-covariance matrix ofθ,p t , t = 1, ..., T provided by a Kalman filter are incorrect due to misspecified serial dependence. The estimated standard errors can be adjusted either by applying the "sandwich" variance estimator, or by using the bootstrap. The bootstrap can additionally adjusts for the non-normality of errors u(t), at the beginning of the epidemic, when the distribution may by closer to a multivariate Poisson distribution than to a normal distribution. The condition for the uniqueness of the solution given in Proposition 2 is an identification condition, which is discussed in detail in the next section. In this section we discuss the (asymptotic) identification corresponding to the objective function with ||p|| 2 = p p given in Section 3. For a homogeneous Markov process with θ = P , this objective function has a simple form, as under linear constraints it is quadratic with respect to the sequence p(t). This allows us for an optimisation in two steps: first with respect to the p(t)'s, and next, with respect to θ after concentrating. This is the approach used below for identification 6 . Next, the analysis is extended to the SIR model to observe the outcomes of a path dependent transmission effect. By taking into account the fact that probabilities sum up to one, we can compare the number of moment conditions equal to (J − 1)(T − 1) with the number of parameters of interest (J −K −1)T + dim θ. Therefore, the order condition is KT −(J −1) ≥ dim θ. It is satisfied iff the number of days T is sufficiently large. However, in a non-linear framework the order condition is insufficient for identification, in general. Let us now consider the rank condition, which is a condition of local identification. For ease of exposition, we first consider the example of a homogenous Markov model with 3 states: J = 3, which is the number of states in a SIR model [see, Section 4.3]. The parameter θ includes the elements of the transition matrix P , which has 6 independent components, given that each row of P sums up to 1. We assume that the observed marginal probabilities are p 3 (t), t = 1, .., T . Thus, we have partial observability. From the Bayes' theorem, it follows that p(t) = P p(t − 1), t = 2, ..., T, (4.5) leading to 2(T − 1) independent moment restrictions that are the estimating equations: p 2 (t) = p 12 p 1 (t − 1) + p 22 p 2 (t − 1) + p 32 p 3 (t − 1), or equivalently, To discuss identification, we search for the solutions in θ = P and p(t), t = 1, ..., T of system (4.3) written for t = 2, ..., T . We have the following result: T ≥ 6, we have that: i) Parameter P is not identifiable, with an under-identification order equal to 3. ii) There exist 3 functions of P that are identifiable. These functions are independent of T . iii) These functions are over-identified with an over-identification order equal to T − 5. Proof: The proof is based on a concentration with respect to the values of p 2 (t). From the second equation of system (4.3), we see that p 2 (t − 1) is a linear affine function of p 3 (t), p 3 (t − 1), with coefficients that depend on P . These linear affine expressions can be substituted into the first equation of system (4.3) to show that the observed sequence p 3 (t) satisfies a linear affine recursion of order 2: with coefficients that depend on P . The results follow since: i) the functions a(P ), b(P ), c(P ) are identifiable; ii) the degree of under-identification of P is: 6-3=3; iii) the degree of over-identification of the identifiable parameters is: T −2−3 = T −5. Q.E.D. Appendix 3 provides the expressions of functions a(P ), b(P ), c(P ) and points out that Proposition 3 holds, except for conditions that are (Lebesgue) negligible. In particular, identification requires that observations p 3 (t) correspond to a nonstationary episode as shown in the remark below. Remark 1: Let π denote the stationary probability solution of the Markov chain, defined by: If the observed p 3 (t) = π 3 were associated to a stationary episode, the sole identifiable function of parameters would be π 3 (P ) and the under-identification degree would be equal to 6-1=5. Therefore, by observing the process during a nonstationary episode, we gain 2 identification degrees. Remark 2: If the Markov structure is recursive, that is, if matrix P is upper triangular, the under-identification degree becomes 3-3=0, and the parameter is generically identifiable. Proposition 3 shows that we can expect to identify the parameter of interest if we either consider a) a homogeneous Markov and constrain the parameters, as illustrated in Remark 2 by an example of the recursive system, or b) a non-homogeneous Markov discussed in the next subsection. Remark 3: The rank condition can be derived in the general case of any number of states J and any type of partial observability of A. The relation between the observations A t (for N large) and the parameters of interest P, p(t), t = 1, ..., T is given by: The second equation can be solved for p(t) as a function of P and p(1), as p(t) = (P ) t−1 p(1). Next, this expression of p(t) can be substituted into the measurement equation to get: Next, we need to find the Jacobian of the transformation associating A 1 , ..., A T to P, p(1). This Jacobian can be obtained by considering the impact of small shocks δP and δp(1) to P and p(1) on A t . By differentiating equation (4.8), we get a linear system in δP and δp (1): System (4.9) can be rewritten in terms of [vec (δP ), vec δp(1)] as: Journal Pre-proof and the rank of Jacobian J can be compared with the parameter dimension (taking into account the unit mass restrictions). In applications, the rank condition has to be checked for each specific model of interest, as shown above for J = 3. Let us now consider an epidemiological model with J = 3 states to facilitate the comparison with the example in Section 4.2. The states of the SID model are: 1=S for susceptible, 2=I for infectious (individuals stay infectious, even if they recover), 3=D for deceased. The rows of the transition matrix are the following: ] is the logistic function, i.e. the inverse of the logit function. We obtain a triangular transition matrix with state D as an absorbing state. The contagion effect is characterized by parameter a 2 and follows a nonlinear logistic function. We also expect that mortality rate p 23 is strictly larger than mortality rate p 13 . There are 4 independent parameters in θ = [a 1 , a 2 , p 13 , p 23 ] . Proposition 4: The SID model with observed p 3 (t) given above is generically identifiable. Parameter θ is over-identified with an over-identification order equal to 5. Proof: The proof is similar to the proof of Proposition 3. The two independent moment conditions are: From the second equation, it follows that p 2 (t − 1) is a linear affine function of p 3 (t) and p 3 (t − 1). Next by substituting into the first equation, we find that the observed p 3 (t) satisfies a nonlinear recursive equation of order 2 of the type: If T is sufficiently large, this nonlinear observed dynamics allows us to identify 9 nonlinear functions of parameter θ. Thus, parameter θ is identifiable with an over-identification order equal to 5. Q.E.D. Remark 2 suggested earlier that the triangular form of the transition matrix alone would facilitate the identification. However, the order of over-identification reveals the additional role of the contagion effect. The nonlinear dynamics induced by the logistic transformation also facilitates identification. Remark 4: As in the case of a homogeneous Markov process, it is theoretically possible to compute the Jacobian associating the observed aggregates A t to the underlying parameters θ, p(t), t = 1, ..., T . The condition on the rank of the Jacobian is difficult to interpret in epidemiological terms, except for specific models, such as the SID model given above. This section illustrates the estimation approach and its performance in an epidemiological model. It is intended to recover the rate of infected undetected individuals, who are often asymptomatic. We consider a model with 5 states: 1=S, 2=IU, 3=ID, 4=R ,5=D, and the following rows of the transition matrix: where the π 1jt , j = 1, 2, 3 sum up to 1, and are proportional to: J o u r n a l P r e -p r o o f Journal Pre-proof McFadden (1984) ]. The transmission parameters b 1 , c 1 , b 2 , c 2 are non-negative and allow for different impacts of p 2 (t − 1) and p 3 (t − 1), as the detected individuals are expected to be self-isolated more often. There is no contagion effect from the recovered R, who are assumed no longer infectious 7 . The structure of zeros in the transition matrix indicates that one cannot recover without being infected, one cannot be infected twice 8 and death is considered as an absorbing state. This is a parametric model with 6+7=13 parameters, i.e. the 6 parameters a l , b l , c l , l = 1, 2 and 7 independent transition probabilities. Among the 5 series of frequencies f j (t), j = 1, .., 5 that sum up to 1 at each date, f 3 (t) and f 5 (t) of infected detected and of deceased, respectively, are assumed to be observed. The frequencies f 2 (t) and f 4 (t) are unobserved and will be considered as additional quantities of interest to be estimated jointly. They are crucial for a model-based inference on counts of infected undetected and of recovered immunized individuals. As illustrated in Section 4.3, the triangular form of the transition matrix and the nonlinear doubly logistic contagion dynamic will provide generic identification. The above model can be used for simulation of the Covid-19 transmission for given values of parameter θ and initial value p(1). These values are set as follows: The daily mortality rates are: We assume that there are about 3 times more transitions to IU than to ID,i.e. The counts of individuals in the two Infected states IU and ID are flows, as they are observed between the times of entry in, and exit from the state of infection. Moreover, the probability of exiting after 20 days is very close to 1. We usually expect a "phase" transition effect: For small t, these counts increase quickly as new infected individuals are cumulated without a sufficiently high number of exits to compensate for the arrivals. This explains an increase of the curves at the beginning of the period. After that initial period, the counts of exits tend to grow and offset the new arrivals so that the curves tend to flatten. More precisely, they continue to increase, due to the disease transmission effect, but at a very low rate. This is the so-called flattening of the curve. This theoretical evolution depends on the choice of parameter values, especially the transmission parameters. Given the selected parameter values that allow for exogenous sources of infection, the initial convex pattern in the counts of infected is not visible. Only the concave part of the curve, up to its flat part, is observed. One can perform similar dynamic sensitivity analysis for other credible scenarios. parameters. A self-isolation measure introduced at some point would have changed subsequent evolution. There is first a tendency to reach a flat part on the curve without self-isolation, and then to reach a lower flat part on the curve with self-isolation measures. Therefore, over a longer period, the first flat part can appear as a smoothed peak. If self-isolation measures are lifted afterwards, a second peak of infections is expected, and so on, resulting in a sequence of stop and go [Ferguson et al.(2020) , Gourieroux, Jasiak (2020)]. We focus on the second period which is sufficiently stable for the estimation purpose. The fully observed states are the states ID and D. State ID is assumed equivalent to hospitalization, as the counts of ("confirmed") detected, which are publicly available, are measured with error and are not reliable. This is due to the counts of detected individuals being derived from the PCR test results, while not all tests results may have been recorded, some people could have been tested multiple times, inflating the counts, or people might have not been tested at random, or without an adequate exogenous stratification, which creates a selectivity bias 13 14 . In contrast, the hospitalization data are more reliable and regularly updated. State D is assumed observed through total death counts. These include deaths from Covid-19, which are reported on-line as D/H, i.e. death after hospitalization, and are known to underestimate the true number of deaths due to the coronavirus, as they do not include all deaths from Covid-19 at home, or in the long-term health care institutions. The series to be estimated are the theoretical proportions of infected undetected IU and recovered R. We use the available series of ("confirmed") detected and of recovered after hospitalization, for comparison with the estimates. More specifically, we use the French data on the total daily number of deaths from the French National Statistical Institute INSEE (2020) and the daily data on coronavirus pandemic from Sante Publique France (2020) reported at https://dashboard.covid19.data.gouv.fr/ and https://www.linternaute.com/actualite/guide-vie-quotidienne/2489651-covid-19-en-franceles-dernieres-statistiques-au-06-avril-2020/, available on April 06 15 . The daily evolutions of total counts of hospitalized, detected, recovered and deceased individuals reported by these sources on April 6 are displayed in Figure 3 . Note that the data used in this study can differ from the data currently reported, due to updating. In particular, the daily data on overall death counts in France have been since updated and adjusted for individuals deceased at home or in long term health care facilities. For example, the new records report 2713 deaths on April 6, 2020, as compared to the initially reported number of 2401 used in this study. if the reproductive number R 0 is larger than 1, and it does not, otherwise. The specification outlined in Section 5.1 differs from the standard SIR in terms of the expressions of transmission functions π 12,t , π 13,t . They are equal to exp(a 1 ), exp(a 2 ), respectively, if p 2 (t − 1) = p 3 (t − 1) = 0, whereas in the standard SIR, they are equal to 0. The estimated non-zero values of exp(â 1 ), exp(â 2 ) reflect the transmission due to individual travelling between countries and regions. The projected results need to be interpreted with caution, due to the uncertainty on parameter estimates [see, Table 2 ]. Another pessimistic outcome is that without any social distancing measures, medical treatment for Covid-19, or a vaccine, it takes about 25 years for the marginal probabilities of IU and ID to decline to 0. This paper is intended to provide a solution for incomplete counts of infected and undetected individuals and of recovered individuals. These unknown quantities can be estimated jointly with the parameters of a compartmental epidemiological model. This approach is illustrated in an estimation involving French count data on Covid-19 infections [see also Brown et al. (2020) for an application to North Carolina]. Our methodology 16 The estimation performed on April 06 could not take into consideration the retrospectively updated total death counts. This could explain, at least to some extent, the observed bias. According to the updated sources, the evolution of deaths was more explosive at the beginning, i.e. close to March 16, and its inflection changed earlier too. required daily data on the total counts of deaths, comprising the deaths due to Covid-19. These data are available in France and other European countries [see, the website Euromomo], but may be publicly unavailable in other countries, such as Canada. The results derived for one country (France) cannot be extrapolated directly to another country or state, because of differences in age structure and comorbidity. More specifically, our results cannot be directly compared with other studies of undocumented infections in the US [see, Hortacsu, Liu, Schwieg (2020) ] and China [Li et al. (2020) ]. The comparisons are difficult, as each study employs different models, aggregate data and estimation methods. For example , Li et. al. (2020) Allen (1994) ]. This difficulty, due to the sensitivity of nonlinear dynamics with respect to the size of timestep, is out of the scope of this paper. Various extensions of the model examined in this paper can be considered: i) As mentioned earlier, the model is a special case of a nonlinear pseudo state space model, with states p(t), deterministic state equations (3.1), and measurement equations: iii) Other specifications of the propagation functions π t can also be considered and compared [see Wu et al. (2020) ]. The treatment of missing data can likely be improved by introducing additional explanatory variables that are expected to impact the virus trnsmission. This approach is followed in Hortacsu et al. (2020) who use hospitalization data from various regions and interregional transportation data to forecast infection rates. where P (t − 1) denotes the transition matrix from date t-1 to date t. By the iterated expectations theorem, we get: Let us now consider the covariance: (by the iterated expectation and using E(Z t ) = p(t)) (by taking into account the 0-1 components of Z) This is the expression of the autocovariance as a function of the p(t)'s and model parameters. Under Assumptions A.1. and after a normalization by 1/N we obtain the autocovariance of the frequencies f (t), t = 1, ..., T and of the measurement equation error p 2 (t) = p 12 + (p 22 − p 12 )p 2 (t − 1) + (p 32 − p 12 )p 3 (t − 1) and substitute into this equation the expression of p 2 (t) derived in part i). We get: It follows that: a(P ) = p 12 (p 23 − p 13 ) + p 13 (1 − p 22 + p 12 ), b(P ) = p 22 − p 12 + p 33 − p 13 , c(P ) = (p 22 − p 12 )(p 13 − p 33 ) + (p 23 − p 13 )(p 32 − p 12 ). To get a recursive equation of order 2, we need the second condition: condition 2: c(P ) = 0 To identify functions a, b, c from the observed p 3 (t), we need: condition 3: The matrix 3×(T −2) with columns (1, ..., 1) , (p 3 (T −1), p 3 (T −2), .., p 3 (2)) and (p 3 (T − 2), p 3 (T − 3), ..., p 3 (1)) is of full column rank. This implies, in particular, the order condition: T ≥ 5 in Proposition 3. The following condition 4 is needed for computing the exact under-identification order of P from functions a, b, c. condition 4: By taking into account the unit mass restrictions on the rows of P ,the Jacobian of (a, b, c) has rank 3. Note that condition 4 implies condition 2. Some Discrete-Time SI, SIR and SIS Epidemic Models A Simple Planning Problem for COVID-19 SIR Model with Time Dependent Infectivity Parameter: Approximating the Epidemic Attractor and the Importance of the Initial Phase Mathematical Models in Population Biology and Epidemiology Time Series Analysis via Mechanistic Models Estimating Undetected COVID-19 Infections. The Case of North Carolina Voluntary and Mandatory Social Distancy: Evidence on COVID 19 Exposure Rates from Chinese and Selected Countries Multinomial Goodness-of-Fit Tests Capturing the Time Varying Drivers of an Epidemic Using Stochastic Dynamical Systems Robust Extended Kalman Filtering Estimating the Number of Infections and the Impact of Non-Pharmaceutical Interventions on Covid-19 in 11 European Countries Competition Between Growths Governed by Bernoulli Percolation Estimating Equations in the Presence of a Nuisance Parameter The Number of Individuals in a Cascade Process Analysis of Virus Transmission: A Transition Model Representaton of Stochastic Epidemiological Models SIR Model with Stochastic Transmission Percolation Processes: Lower Bounds for the Critical Probability The Mathematics of Infectious Diseases Estimating the Fraction of Unreported Infections in Epidemics with a Known Epicenter: An Application to COVID-19 Nombre de deces quotidiens par departement A New Extension of the Kalman Filter to Nonlinear Systems,11th Int.Symp. on Aerospace/Defence,Sensing,Simulation and Controls A Contribution to the Mathematical Theory of Epidemics Directions in Mathematical Systems:Theory and Optimization Substantial Undocumented Infection Facilitates the Rapid Dissemination of Novel Coronavirus(SARs-CoV2) Estimating the COVID-19 Infection Rate:Anatomy of an Inference Problem Econometric Analysis of Qualitative Response Models Estimation of Time Varying Markov Processes with Aggregate Data Applications of Mathematics to Medical Problems Information Recovery in a Dynamic Statistical Markov Model, Econometrics Estimation of the Asymptomatic Ratio of Novel Coronavirus Infection (COVID-19) Donnees Hospitalieres Relatives a l'Epidemie Covid-19 The Extended Kalman Filter as a Local Asymptotic Observer Estimation of the Transmission Risk of the 2019-nCoV and its Implication for Public Health Interventions Estimates of the Severity of Coronavirus Disease 2019;A Model Based Analysis An Introduction to Infectious Disease Modelling The Unscented Kalman Filter for Nonlinear Estimation Generalized Logistic Growth Mod The asymptotic expansions are easily derived, given that the optimization in Proposition 2 is deterministic. Therefore, estimatorsp(1),p(2), ...,p(T ),θ are deterministic functions of observations t = Af (t), t = 1, .., T . If the transition matrix is twice continuously differentiable with respect to p(t − 1) and θ in a neighbourhood of the true values, these deterministic functions are continuously differentiable. Then, by using the asymptotic normality of f (t)'s (Proposition 1), we can apply the delta method to deduce the 1/ √ N rate of convergence of the estimators and their asymptotic variance-covariance matrix from the one of the f (t)'s (see Appendix 1).When the number of observation dates and of missing counts is too large, the use of the delta method can be numerically cumbersome. It can be replaced by a bootstrap method (for which the regularity conditions of validity are satisfied in our framework), or by the approximated standard errors provided by an EKF, or UKF algorithm, after adjusting for the misspecification of the autocovariances of the measurement equation errors u(t).Appendix 3 This Appendix derives the equations used in the proof of Proposition 3. It provides the closed form expressions of functions a(P ), b(P ), c(P ), and outlines conditions 1 to 4for the validity of Proposition 3.i) Let us first solve the second equation of system (4.3). We get: