Forecasting Elections from Partial Information Using a Bayesian Model for a Multinomial Sequence of Data

Voter models have been of interest for a long time. With the strengthening of democracy and the role of media, statistical models for the prediction of elections have been on the rise around the world. Elections are studied at different stages – starting with pre-poll predictions, going on to post-poll predictions, and also calling the election before the final result is declared as the counting of ballots is in progress. The current work is focuses on the last context, where the counting of results is publicly disclosed in batches (or rounds). Stemmed from public interest, this often leads to competition among the media houses “to call the election”, i.e. to predict the outcome correctly and as early as possible.

Using the exit poll data has been a popular approach in this context. It however has its own pitfalls, as has been shown in many papers (see, e.g. Stiers and Dassonneville [2018]). Because of the inadequacy of the exit poll data, Associated Press (AP) designed its methodology for both survey and forecasting (Slodysko [2020]). Interestingly, they still do not call a closely contested race unless the trailing candidate cannot mathematically secure a victory. Another common approach is to extract relevant information from social media, such as Twitter or Google Trends, and to use that to predict the winner of an election. O’Connor et al. [2010] and Tumasjan et al. [2010] are two notable works in this regard. A major criticism behind this approach is that intentional bias in the social media posts often lead to error-prone conclusions in the forecasting problems (Anuta et al. [2017]). Further, Haq et al. [2020] provides a great review of the existing literature on how to call an election and one can see that most of the methods are either based on data from social media or are ad-hoc and devoid of sound statistical principles. To that end, it would be of paramount importance to develop a good forecasting methodology based only on the available information of votes secured by the candidates in the completed rounds of counting. This can be treated as a collection of multinomial data where the number of categories is equal to the number of candidates and the random fluctuations in cell probabilities point to the randomness of the vote shares in different rounds.

Multinomial distribution is the most common way to model categorical data. In this work, we study the behaviour of a collection of multinomial distributions with random fluctuations in cell (or ‘category’, as we interchangeably use the two terminology) probabilities. We consider situations where this multinomial data is observed in batches, with the number of trials in each batch being possibly different. The cell probabilities for every multinomial random variable in the data are assumed to be additive combinations of a fixed component which is constant across batches, and a random perturbation. These perturbations for different batches are taken to be independent. Before delving deeper into the model and related results, we present an interesting application where the proposed setup can be utilized effectively.

In this section, we consider a few toy election scenarios and evaluate the effectiveness of our proposed methodology across various scenarios of elections. We generate $K$ batches of election data as multinomial random draws following eq. 2.2 for three different simulated election outcomes (SEO). Each election data is assumed to have $C$ candidates. The expected values of $n_{j}$ (the number of votes cast in the $j^{th}$ batch ) are considered to be the same for all $j$. Let $n$ denote the common expected value for different batches. We consider different value of $C$ (3 or 5), $K$ (25 or 50), $n$ (100 or 1000 or 5000 or 50000) to understand the robustness of the method.

The SEOs we consider differ in the way the $\bm{\varepsilon}_{j}$’s are generated. First, we take the simplest situation of $\varepsilon_{cj}=0$ for all $c,j$, which corresponds to an iid collection of multinomial random variables. Second, $\bm{\varepsilon}_{j}$ is simulated from a truncated multivariate normal distribution with mean $\bm{0}$ and dispersion matrix $n_{j}^{-0.5}\mathbf{I}$. Third, we impose dependence across different coordinates of $\bm{\varepsilon}_{j}$, and simulate it from a truncated multivariate normal distribution with mean $\bm{0}$ and dispersion matrix $n_{j}^{-0.5}A$, where $A$ is a randomly generated positive definite matrix with finite norm. Observe that all the three SEOs satisfy 1. For every combination of $K,C,n$, each experiment is repeated many times and we compute the average performance over these repetitions.

In this simulation study, we particularly explore two aspects in detail. First, we find out the accuracy of our methodology in terms of predicting the candidate with maximum votes at the end. Here, a decision is made when the winner is predicted with at least 99.5% probability and when the predicted difference in votes for the top two candidates is more than 5% of the remaining votes to be counted. Second, we focus on two specific cases which are similar to our real-life applications, and examine how well we can predict various properties of $h(\bm{Y}_{l})$, for different choices of $h(\cdot)$. One of these cases has $C=5,K=50$ and we look for different options for $n$. In the other case we have $C=3,K=25,n=5000$.

In case of the first problem, for every iteration, the fixed values of $\bm{p}$ are obtained from a Dirichlet distribution with equal parameters to ensure a more general view of the performance of the proposed method. Detailed results are reported in Appendix A, in Tables A1, A2, and A3. We note that the results are not too sensitive with respect to the choices of $C$ and $K$, but the accuracy increases with $n$. Below, in Table 1, we look at the results for $C=5$, $K=50$ for the three different SEOs. It is evident that under the iid assumption (SEO1), the method makes a correct call almost all the times. When the components of $\bm{\varepsilon}_{j}$ are independent random variables (SEO2), a correct call is made more than 90% of the times for large enough samples ($n>5000$). For the third SEO, where the components of $\bm{\varepsilon}_{j}$ are simulated from a multivariate distribution with a general positive definite matrix as the covariance matrix, the accuracy of making a correct call drops to 70% for $n=50000$. Another interesting observation is that the accuracy depends heavily on the average final margin between the top two categories. In other words, if the underlying probabilities of the top two candidates are close (which directly relates to a lower margin of difference between the final counts), then the method achieves lower accuracy, and vice-versa. We also observe that in those cases, the method records no call (i.e. a decision cannot be made with the prescribed rules) more often than incorrect calls.

Moving on to the particular choice of $C=3,K=25,n=5000$, we aim to find out the effectiveness of the proposed approach in identifying the leading candidates for varying degrees of difference between the true vote probabilities of the topmost two candidates. To that end, for $\bm{p}=(p_{1},p_{2},p_{3})$, without loss of generality, we assume $p_{1}>p_{2}>p_{3}$. Now, data are generated for the above three simulated elections by fixing $p_{1}-p_{2}=\delta$ where $\delta\in\{0.01,0.05,0.1,0.25\}$. For every choice of $\delta$ and for every SEO, we repeat the experiments and find out the mean accuracy in predicting the top candidate. Additionally, we also find out the average number of observations used by the method for making the prediction. These results are displayed in Table 2.

It is evident that for SEO1, the method provides great prediction. Even when there is only 1% difference in the probabilities of the top two candidates, the method predicts the leading candidate correctly 96% of the times. For larger than 1% difference, it never fails to predict the leading candidate correctly. For SEO2, our method can predict the topmost candidate with more than 75% accuracy whenever the probabilities of the top two candidates differ by at least 5%. This accuracy reaches the value of nearly 90% for $\delta=0.1$ and is 100% if the difference in those two cell probabilities is 0.25. For the third SEO, the prediction accuracy is about 70% for $\delta=0.1$. The accuracy improves steadily as $\delta$ increases. One can say that the method will be able to predict the winner with high level of accuracy whenever there is a considerable difference between the $p_{i}$ values for the top two candidates. If that difference is minute, which corresponds to a very closely contested election, the accuracy will drop. Additionally, we observe that in the most extreme cases, about 20% of the times, the method never declares a winner with desired certainty, which is in line with what we should expect.

We now focus on the second specific case of $C=5,K=50,n=50000$. For different SEOs, the values of $\bm{p}$ are generated randomly from a Dirichlet distribution, and we evaluate how well our method can predict the cumulative proportions of counts for the five categories. Experiments are repeated many times and we compute the overall root mean squared error (RMSE) for every SEO, based on the predictions made at different stages. Refer to Table 3 for these results. Once again, for SEO1, the method achieves great accuracy very early. With only 15 rounds of data (approximately 30% of the total observations), the predictions are precise. At a similar situation, for SEO2, the RMSE is less than 1%, and it decreases steadily to fall below 0.5% by round 35. Finally, for the third SEO, RMSE of less than 1% is recorded after 35 rounds as well. These results show that our method can accurately predict the overall cumulative proportions of counts with about 35 rounds of data. Translating this to the sales forecasting problems, we hypothesize that the proposed method will be able to predict the annual market shares of different categories correctly at least 3 to 4 months before the year-end.

We examine the effectiveness of our proposed methodology in the election calling context, with the real data from the Bihar Legislative Assembly Election held in 2020, across 243 constituencies (seats). Indian elections have multi-party system, with a few coalition of parties formed before (sometimes after) the election, playing the major role in most contests. While these coalitions are at times temporary in nature and there is no legal binding or legitimacy, typically the alliance winning the majority (at least 50%) of the constituencies comes to power. In the Bihar 2020 election, the contests in most constituencies were largely limited to two such alliances – the National Democratic Alliance (NDA), and the Mahagathbandhan (MGB). Although there were a couple of other alliances, they did not impact the election outcomes and for the purpose of this analysis, we consider them in “others” category. For each constituency, the counting of votes took place in several rounds, and we want to examine how these round-wise numbers can be useful in predicting the final winner ahead of time.

The result of the election is downloaded from the official website of Election Commision of India [2020]. A brief summary of the data is provided in Table 4. All summary statistics in this table are calculated over the 243 constituencies in Bihar (our analysis ignores the postal votes, which are few in numbers and do not affect the results). Note that the number of votes cast in the constituencies vary between 119159 and 225767, with the average being 172480. The votes are counted in 25 to 51 rounds across the state, and the mode is found out to be 32. While in most rounds about little over 5000 votes are counted, in some cases, especially the last few rounds, show fairly less numbers. The last three rows in the table provide the minimum, maximum and average number of votes counted per round for all the constituencies. We also point out that the final margin between winner and runner-up parties in these constituencies show a wide range. The lowest is recorded in Hilsa (13 votes) whereas the maximum is recorded in Balrampur (53078 votes).

As mentioned above, in the following analysis, we focus on the two dominant categories (essentially the two main alliances in every constituency) and the rest of the alliances or the parties are collated into the third category. As our main objective is to predict the winner and the margin of victory, it is sensible to consider this three-category-setup. We also make the logical assumption that each vote is cast independently, thereby ensuring that in each round of counting, we have an independent multinomial distribution. However, there is no information on how the counting happens sequentially in each round, and therefore we assume the individual probabilities of the three different categories in different counting rounds to be random variables satisfying 1. Consequently, it is an exciting application of the proposed modeling framework.

Recall the hierarchical prior distributions for the parameters in the model (eq. 2.11). For $\bm{\alpha}$, we use the past election’s data and use the proportion of votes received by the corresponding alliances (scaled, if needed, to have $\sum\alpha_{i}=1$). For $\Psi$ and $\Psi_{p}$, we use identity matrices of appropriate order. Both $\nu$ and $\nu_{p}$ are taken to be 5. We point out that the results are not too sensitive to these choices. Next, using these priors, we implement our method and estimate the probability of winning for every alliance after the counting of votes in each round. Corresponding margin of victory is also estimated in the process. Based on these values, we propose the following rules to call a race. First, we allow at least 50% votes to be counted in order to make a prediction. It is otherwise considered to be too early to call. On the other hand, akin to the procedure laid out by Lapinski et al. [2020], our decision is to call the race in a particular constituency when we are at least 99.5% confident of the winner and when the predicted winning margin is at least 5% of the remaining votes. Unless these conditions are met, it is termed as too close to call.

We start with a summary of the results obtained after fitting the model to the election data. Complete results are provided in LABEL:tab:bihar-fullresults in Appendix A. In 227 out of the 243 constituencies (approximately 93.4%), the proposed method calls the race correctly. Our method considers that it was always too close to call for Bakhri and Barbigha where the eventual win margin was 439 and 238 respectively. For 14 other constituencies (approximately 5.7%), the method calls the race incorrectly, as the winning candidate was significantly behind at the time of calling in those constituencies.

We next take a detailed look at the prediction patterns for all constituencies, in different aspects, through Figure 2. In both panels of the Figure we color the correct predictions by green, the incorrect predictions by red, and no calls by blue. In the left panel of the Figure we plot the remaining vote percentages at the time of calling along the color gradient. In the right panel we plot the final margins along the color gradient. It is interesting to observe that in the right panel all red points, which correspond to incorrect calls, appear with lighter shades. Thus, our method performs perfectly for all constituencies but one where the final margins of victories are above 10000. From the left panel it is also evident that the proposed method calls the races very early for all such constituencies. Only exception to this is the constituency Baisi, which is discussed in more detail below. We can also observe that seats painted in lighter shades in left panel are also in lighter shades on right panel, implying that the closely contested constituencies require more votes to be counted before the race can be called.

To further explore the above point, the accuracy and the average percentage of votes counted before making a call, corresponding to the final margins of victories, are displayed in Table 5. There are 52 constituencies which observed very closely fought election (final margin is less than 5000 votes). In 41 out of them, our method correctly predicts the winner and in 2 of them, it cannot make a call. On average, around 70% votes are counted for all these constituencies before we can make a call. This, interestingly, drops drastically for the other constituencies. On average, only about 60% votes are counted before we make a call in the constituencies where eventual margin is between 5000 and 10000. In 4 out of 34 such constituencies though, our method predicts inaccurately. For the constituencies with higher eventual margins (greater than 10000), the prediction turns out to be correct in 156 out of 157 constituencies. In these cases, around only about 50 to 55 percentage of votes are counted on average before making a call. Finally, the last row of the table corroborates the earlier observation that the size of the constituency does not have considerable effect on the prediction accuracy.

Not only the winner of an election, but the proposed approach also predicts the final margin of victory, along with a prediction interval, for the winner. In Figure 3, we present the true margin and the predicted win margin for all the 227 constituencies where correct calls are made. Corresponding prediction intervals are also shown in the same plot. It can be observed that the predicted margins closely resemble the true margins.

Overall, it is evident that the method performs really well in terms of calling an election race based on only the counting of votes in different rounds. It is imperative to point out that better accuracy can be achieved with more detailed information about different constituencies. Especially, information about the sequence in which the counting is carried out, ethnic composition or socio-economic status of people in different areas are necessary to build a more accurate procedure. Our proposed method in this paper works with only the multinomial assumption and achieves more than 93% accuracy in calling the race with the requirement of at least 50% votes being counted. To that end, in order to gain more insight about the accuracy of the method, we rerun the experiment with different values for the minimum requirement of counting. The number of constituencies with correct call, incorrect call and no call are displayed against the minimum requirement in Table 6. We see that the method achieves 90% accuracy (215 out of 243) even when only 30% votes are counted. For minimum counting of 60% or more votes, the results do not improve substantially, except for the fact that no call can be made in more number of constituencies.

Next, we look into the robustness of the method, by focusing on four particular constituencies in more detail. These are Hilsa (minimum final margin), Balrampur (maximum final margin), Baisi (final margin more than 10000, but our proposed methodology makes a wrong call with nearly 55% votes counted) and Barbigha (final margin of 238, the proposed method never calls the race in anyone’s favour). Now, for each constituency, we create synthetic data by randomly switching the order of the individual rounds of counting and implement our method on the resulting data. This exercise helps us to understand whether the predictions would have been different if the data were available in a different order. In particular, it points towards objectively judging to what extent the success or failure of the proposed methodology in this real data can be attributed to luck. Note that the proposed modeling framework assumes independence across different rounds, and hence the method should perform similarly if the votes are actually counted in random order. In Table 7 below, the results of these experiments are presented. We see that for Balrampur, even with randomly permuted rounds, the method calls the race correctly every time. For Baisi, the method yields correct forecast 99.97% of times, although in the original order of the data it makes an incorrect call. Meanwhile, Hilsa and Barbigha observed very closely contested elections. For these two constituencies, in the experiment with permuted data, the method records too close to call decisions more commonly, 72.55% and 70.53% respectively. These results suggest that quite possibly the sequence of roundwise counting of votes is not random. Thus, if more information is available about individual rounds from which votes are counted, it would be possible to modify the current method to achieve greater accuracy.

In this section, in an attempt to impress upon the reader that the proposed methodology can be extended to other election forecasting contexts as well, we focus on the statewise results from the 2020 Presidential Election of the USA (data source: MIT Election Data and Science Lab [2020]). It is well known that in American politics, every election relies primarily on a few states, popularly termed as the ‘swing states’, where the verdict can reasonably go to either way. For the other states, the results are guessable with absolute certainty, and are of less interest in any such research. Thus, in this application, our focus is on the eleven swing states of the 2020 Presidential Election – Arizona, Florida, Georgia, Michigan, Minnesota, Nevada, New Hampshire, North Carolina, Pennsylvania, Texas and Wisconsin. Since the final outcome of the election in these cases are extremely likely to be either democratic or republican, with a minimal chance of a third-party-victory, we can consider a three-category-setup. Then, akin to Section 4, it can be easily argue that the setup follows the proposed modeling framework. We also adopt the same setting for prior distributions and the procedure to call the race.

Our objective is to demonstrate that partial data of the counties can provide accurate forecast of the final result in these turbulent states. However, we note that the information of the particular order in which results were declared for different counties are unavailable. That is why we call it a hypothetical application, and we circumvent the issue by taking a permutation-based approach to implement the proposed method. For every state, we randomly assign the orders to the counties and assume that the data are available in that sequence. This permutation procedure is repeated 100 times, and we evaluate the performance of the proposed method at a summary level across these permutations.

Let us start with a brief overview of the number of counties and the number of voters in the swing states (see Table 8). All summary statistics are calculated based on the eleven states. We can see that the number of counties as well as the number of voters vary substantially; although the final margin of victory remains within 7.4% of the total votes. It is in fact less than 1% in three states (Arizona, Georgia and Wisconsin), making them the most closely contested states in this election.

Turn attention to the results of the forecasting procedure on the permuted datasets. In the last three columns of Table 9, we display the proportion of times our method correctly calls the race, finds it too close to call, or incorrectly calls the race based on partial county-level data,

We observe that in New Hampshire, Florida, Nevada and Texas, partial data provide the correct prediction in almost all cases. One may connect this to the fact that the final margin in these states were on the higher side, especially if we look at the values in terms of proportion to the total number of votes. Arizona, on the other hand, was the most closely contested state with a big voter pool. Our approach is still able to predict the correct winner 60% of the times whereas it decides the race to be too close to call on 39% occasions. Georgia and Wisconsin are the two other states with the lowest final margins. There, our algorithm makes wrong judgements in more than 40% situations.

In the other four states, namely Michigan, Minnesota, North Carolina and Pennsylvania, the results are not good. Albeit the final margins are considerably large, the final conclusions are wrong in at least 18% of the permuted datasets. This phenomena brings forward an important limitation of the methodology. In our procedure, the key assumption is that the perturbations, albeit random and different in different counties, have the same asymptotic distribution. In such political applications, this assumption may be violated on certain occasions. To demonstrate it empirically, let us look at the margins of differences across the counties in these states (Figure 4).

From the figure, it is evident that in all these states (most prominently in Michigan, Minnesota, Georgia and Wisconsin), the final verdict relied heavily on only a couple of counties. Quite naturally, a partial dataset with the information of those particular counties is extremely likely to provide contradictory results to a partial dataset without those counties. In North Carolina, on the other hand, several counties observed high margin of difference in favor of the democratic party, while the final ruling went the opposite way. These situations point to a potential deviation from the aforementioned assumption, and we hypothesize that in such circumstances, the method can be updated to include more information about individual counties. It calls for an interesting future research direction and we discuss it in more detail in the next section.

In this work, we consider a hierarchical Bayes model for batches of multinomial data with cell probabilities changing randomly across batches. The goal is to predict various properties for aggregate data correctly and as early as possible. We illustrate the effectiveness of the methodology in poll prediction problems which motivated this study. The performance of the methodology is good, especially considering the limited information. We plan to implement the method in future elections. A potential future direction in this specific context would be to improve the method when more information on the constituencies or the rounds of counting are available. It is likely that additional information on relevant covariates would improve forecast accuracy. It is possible to modify the proposed methodology to incorporate the covariates, but that is a more challenging and interesting problem and requires full treatment. We defer it for a future work.

The proposed model accommodates randomness in the multinomial cell probabilities, depending on the batch. This is pragmatic given the actual data pattern in most practical situations, as otherwise early calls are made which end up being often wrong. This is also intuitively justifiable as, for example, different rounds of votes can have differential probabilities for the candidates, potentially owing to different locations or the changes in people’s behaviour. This paper aims to tackle these situations appropriately.

It is worth mention that a common practice in Bayesian analysis of sequences of multinomial data is to apply Dirichlet priors for the cell probabilities. This approach requires the iid assumption and hence do not work under the modeling framework of this study. In fact, if we assume iid behaviour for the multinomial data across the batches and use that model for the examples discussed above, then the performances are much worse than what we achieve. Further, we want to point out that even without the variance stabilizing transformation, the proposed Bayesian method in conjunction with the results from Proposition 1 can be used in similar problems. The prediction accuracy of that model is comparable to our approach. We employ the transformation for superior prediction performance, albeit marginally in some cases, across all scenarios.

Finally, note that in some contexts, it may be unrealistic to assume that values of the future $n_{j}$’s are known. In that case, the decision can be taken by anticipating them to be at the level of the average of $n_{j}$’s observed so far. It is intuitively clear from Section 2.4, as one can argue that the effectiveness of the methodology will not alter to any appreciable extent as long as the sample sizes are large. It can also be backed up by relevant simulation studies.

The methodology proposed in this work can be applied in other domain as well, most notably in sales forecasting. Retail giants and companies track the sales for different products in regular basis. This is useful from various perspectives. One of the key aspects is to project or assess the overall (annual) sales pattern from the data collated at monthly, weekly or even daily level. It is of great value to a company to infer about its eventual market share for the period (year) as early as possible, hence to validate if its (possibly new) marketing plans are as effective as targeted. The retailer would find these inferences useful from the perspective of managing their supply chain, as well as in their own marketing. Other possible applications of the proposed methodology (possibly with necessary adjustments) include analyzing hospital-visits data, in-game sports prediction etc. Some of these applications can make use of additional covariate information; we propose to take that up in a follow-up article.

• Proposition 1.Following eq. 2.2, $\mathbb{E}(X_{lj}\hskip 1.0pt|\hskip 1.0pt\varepsilon_{lj})=n_{j}(p_{l}+\varepsilon_{lj})$. Then, using 1 and applying the relationship between total and conditional expectation,

  $$\mathbb{E}(X_{lj})=\mathbb{E}\left[\mathbb{E}(X_{lj}\hskip 1.0pt|\hskip 1.0pt\varepsilon_{lj})\right]=n_{j}p_{l}.$$

  (7.1)

  We also know that $\mathrm{Var}(X_{lj}\hskip 1.0pt|\hskip 1.0pt\varepsilon_{lj})=n_{j}(p_{l}+\varepsilon_{lj})(1-p_{l}-\varepsilon_{lj})$. Then, we can use the relationship $\mathrm{Var}(X_{lj})=\mathrm{Var}\left(\mathbb{E}(X_{lj}\hskip 1.0pt|\hskip 1.0pt\varepsilon_{j})\right)+\mathbb{E}\left(\mathrm{Var}(Y_{lj}\hskip 1.0pt|\hskip 1.0pt\varepsilon_{j})\right)$ to get the following:

  $$\mathrm{Var}(X_{lj})=\mathrm{Var}\left(n_{j}(p_{l}+\varepsilon_{lj})\right)+\mathbb{E}\left[n_{j}(p_{l}+\varepsilon_{lj})(1-p_{l}-\varepsilon_{lj})\right]=n_{j}p_{l}(1-p_{l})+n_{j}(n_{j}-1)\mathrm{Var}(\varepsilon_{lj}).$$

  (7.2)

  Finally, note that the covariance of $X_{1j}$ and $X_{2j}$, conditional on $\bm{\varepsilon}_{j}$, is equal to $-n_{j}(p_{1}+\varepsilon_{1j})(p_{2}+\varepsilon_{2j})$. One can use this and apply similar technique as before to complete the proof of part (a).

  Part (b) follows from (a) using the fact that the multinomial data in different batches are independent of each other.

  To prove part (c), we use the multinomial central limit theorem which states that if $(D_{1},D_{2},\ldots,D_{C})\sim\mathrm{Multinomial}(n,\pi_{1},\pi_{2},\ldots,\pi_{C})$, then

  $$\sqrt{n}\left(\begin{pmatrix}D_{1}/n\\ D_{2}/n\\ \vdots\\ D_{C-1}/n\end{pmatrix}-\begin{pmatrix}\pi_{1}\\ \pi_{2}\\ \vdots\\ \pi_{C-1}\end{pmatrix}\right)\overset{\mathcal{L}}{\rightarrow}\mathcal{N}_{2}\left(0,\begin{bmatrix}\pi_{1}(1-\pi_{1})&-\pi_{1}\pi_{2}&\ldots&-\pi_{1}\pi_{C-1}\\ -\pi_{1}\pi_{2}&\pi_{2}(1-\pi_{2})&\ldots&-\pi_{2}\pi_{C-1}\\ \vdots&\vdots&\ddots&\vdots\\ -\pi_{1}\pi_{C-1}&-\pi_{2}\pi_{C-1}&\ldots&\pi_{C-1}(1-\pi_{C-1})\end{bmatrix}\right).$$

  (7.3)

  Let us define $V_{\bm{\varepsilon}_{j}}$ to be the covariance matrix of $(\hat{p}_{1j},\hat{p}_{2j},\ldots,\hat{p}_{(C-1)j})$, conditonal on $\bm{\varepsilon}_{j}$. Then, for the $(c,c^{\prime})^{th}$ element of $V_{\bm{\varepsilon}_{j}}$, we have

  $$(V_{\bm{\varepsilon}_{j}})_{c,c^{\prime}}=\begin{cases}(p_{c}+\varepsilon_{cj})(1-p_{c}-\varepsilon_{cj})&\text{for }c=c^{\prime},\\ -(p_{c}+\varepsilon_{cj})(p_{c^{\prime}}+\varepsilon_{c^{\prime}j})&\text{for }c=c^{\prime}.\end{cases}$$

  Now, under 1 and from eq. 7.3, as $n_{j}\to\infty$,

  $$\sqrt{n_{j}}V_{\bm{\varepsilon}_{j}}^{-1/2}\left(\hat{\bm{p}}_{j}-\tilde{\bm{p}}-\bm{\varepsilon}_{j}\right)\mid\bm{\varepsilon}_{j}\overset{\mathcal{L}}{\rightarrow}\mathcal{N}_{C-1}\left(\bm{0},\mathbf{I}\right).$$

  (7.4)

  Using the assumptions in part (c) of Proposition 1, as $n_{j}\to\infty$, $V_{\bm{\varepsilon}_{j}}$ converges in probability to $\Xi_{p}$. Hence, by an application of Slutsky’s Theorem, eq. 7.4 is equivalent to:

  $$\sqrt{n_{j}}\;\Xi_{p}^{-1/2}\left(\hat{\bm{p}}_{j}-\tilde{\bm{p}}-\bm{\varepsilon}_{j}\right)\mid\bm{\varepsilon}_{j}\overset{\mathcal{L}}{\rightarrow}\mathcal{N}_{C-1}\left(\bm{0},\mathbf{I}\right).$$

  (7.5)

  Let $f_{\bm{\varepsilon}}(\cdot)$ denote the density function of $\sqrt{n_{j}}\bm{\varepsilon}_{j}$, $\phi(\cdot,\Xi_{\bm{\varepsilon}})$ be the density function of $\mathcal{N}_{C-1}(\bm{0},\Xi_{\bm{\varepsilon}})$, and $\Phi(\cdot,\Xi_{\bm{\varepsilon}})$ be the distribution function of $\mathcal{N}_{C-1}(\bm{0},\Xi_{\bm{\varepsilon}})$. Note that the previous convergence in eq. 7.5 is conditional on $\bm{\varepsilon}_{j}$, as $n_{j}\to\infty$. To find the unconditional distribution of $\hat{\bm{p}}_{j}-\tilde{\bm{p}}$, we look at the following, for $z\in\mathbb{R}^{C-1}$:

  $$\mathbb{P}(\sqrt{n_{j}}(\hat{\bm{p}}_{j}-\tilde{\bm{p}})\leqslant z)=\int_{\mathbb{R}^{C-1}}\mathbb{P}(\sqrt{n_{j}}(\hat{\bm{p}}_{j}-\tilde{\bm{p}}-\bm{\varepsilon}_{j})\leqslant z-y\mid\bm{\varepsilon}_{j})f_{\bm{\varepsilon}}(y)\;dy.$$

  (7.6)

  From the assumptions in part (c) of Proposition 1, $\sqrt{n_{j}}\bm{\varepsilon}_{j}$ converges in distribution to $\mathcal{N}_{C-1}(\bm{0},\Xi_{\bm{\varepsilon}})$. Then,

  $$\lim_{n_{j}\to\infty}\mathbb{P}(\sqrt{n_{j}}(\hat{\bm{p}}_{j}-\tilde{\bm{p}}-\bm{\varepsilon}_{j})\leqslant z-y\mid\bm{\varepsilon}_{j})f_{\bm{\varepsilon}}(y)=\Phi(z-y,\Xi_{p})\phi(y,\Xi_{\bm{\varepsilon}}).$$

  (7.7)

  A simple application of dominated convergence theorem then implies, as $n_{j}\to\infty$,

  $$\mathbb{P}(\sqrt{n_{j}}(\hat{\bm{p}}_{j}-\tilde{\bm{p}})\leqslant z)\to\int_{\mathbb{R}^{C-1}}\Phi(z-y,\Xi_{p})\phi(y,\Xi_{\bm{\varepsilon}})\;dy.$$

  (7.8)

  Note that the term on the right hand side of eq. 7.8 is the distribution function of $\mathcal{N}_{C-1}(\bm{0},\Xi_{\bm{p}}+\Xi_{\bm{\varepsilon}})$ at the point $z\in\mathbb{R}^{C-1}$. Hence we arrive at the required result, eq. 2.5. ∎