The coronavirus disease pandemic has shown profound impacts on public health and the economy worldwide. The development of efficient and effective public health interventions to prevent major outbreaks and contain the pandemic relies heavily on a quantitative understanding regarding the spread of the virus, such as the transmission rate and the average incubation period. A commonly used approach in epidemiology is to estimate these quantities of interest using epidemic mathematical models, such as the susceptible-infected recovered (SIR) model, with agent-based simulations which capture complex social networks and global scale into the models (Funk et al., 2009; Heesterbeek et al., 2015; Epstein, 2009 ).
To estimate the parameters of interest, a widely used frequentist approach is to minimize the sum of the squared distances between the observed data and the simulation outputs, which is often referred to as the least squares approach. See, for example, Chowell et al. (2003 Chowell et al. ( , 2004 ; Capaldi et al. (2012) ; Chowell (2017) ; Anastassopoulou et al. (2020) ; Bentout et al. (2020) ; Chen and Qiu (2020) ; Giordano et al. (2020) . This estimation approach is intuitive and easy to compute; however, it is shown in this paper that this method is inefficient, particularly when the mathematical models associated with the simulators are built under certain assumptions or simplifications, which may not hold in reality. These simulators are called imperfect simulators in the computer experiment literature (Kennedy and O'Hagan, 2001; Tuo and Wu, 2015; Plumlee, 2017) . Imperfect simulators are common in epidemiology (Heesterbeek et al., 2015) , and therefore the parameter estimates based on the least squares approach are not efficient.
To improve the estimation efficiency, in this paper we propose a new estimation method for the parameters of interest in the epidemic models. In the computer experiment literature, these unknown parameters associated with the mathematical models are often called calibration parameters, and the process of estimating the parameters such that the model simulations agree with the observed data is called calibration (Kennedy and O'Hagan, 2001; Santner et al., 2018) . Although there are numerous developments on calibration, most of the work focus on continuous outputs while the discussions on non-Gaussian outputs, such as count data which are often observed in epidemiology, are scarce (Sung et al., 2020a; Grosskopf et al., 2020) . In this paper, we propose a new estimation method for non-Gaussian outputs, particularly for count data for our applications in epidemiology, which minimizes the L 2 projection of the discrepancy between the true mean process and the simulation outputs. It can be shown that the proposed estimator is asymptotically consistent, and provides a smaller asymptotic variance than the least squares estimator. Furthermore, it can be shown that the proposed estimator achieves the semiparametric efficiency, even when the model simulations cannot match the reality due to certain assumptions or simplifications.
It is worth noting that there are extensive studies and applications of calibration by Bayesian procedures (Diekmann et al., 2013; Farah et al., 2014; Wang et al., 2020; Wu et al., 2020) . However, without taking the model imperfection into account in the conventional Bayesian framework, the theoretical justification for the parameter estimation with imperfect simulators are not fully developed. On the other hand, Bayesian calibration of Kennedy and O'Hagan (2001) takes into account the model imperfection through Gaussian process modeling, but it suffers from the unidentifiability issue when the parameter estimation is of interest (Bayarri et al., 2007; Han et al., 2009; Hodges and Riech, 2010; Paciorek, 2010; Gramacy et al., 2015) . Furthermore, most of the existing developments are based on continuous outputs with a Gaussian assumption, which is not valid for the count data in the epidemic models in our applications. Recent studies on addressing the unidentifiability issue can be found in Plumlee (2017) and Tuo (2019) .
The remainder of the paper is organized as follows. In Section 2, the new estimation method based on L 2 projection is introduced. Theoretical properties of the proposed estimator are developed in Section 3. In Section 4, numerical studies are conducted to demonstrate the finite sample performance of the proposed estimator and the empirical comparison with the least squares estimator. In section 5, the estimation method is applied to the study of COVID-19 based on an SEIR (Susceptible-Exposed-Infectious-Recovered) compartmental model with both deterministic and stochastic simulations. Discussions and concluding remarks are given in Section 6. Computational details for the estimation are given in Appendix, and the mathematical proofs and the R (R Core Team, 2018) code for implementation are provided in Supplementary Material.
2 Estimation for Compartmental Models in Epidemiology 2.1 Least squares estimator Compartmental models are widely used mathematical models for the study of infectious diseases, in which the population is assigned to compartments with labels such as S (Susceptible), I (Infectious), and R (Recovered and immune) (Diekmann et al., 2013) . For example, the most basic compartment model in epidemiology is the SIR model, in which each living individual is assumed to be in one of the above three compartments at any given time. The model specifies the rates at which individuals change their compartment, which include the progresses from S to I when infected, and from I to R upon recovery.
Various extensions to SIR have been extensively developed in the literature of epidemiology (Heesterbeek et al., 2015) .
Let f (x, θ) denote the number of people in the infectious compartment at time x ∈ Ω ⊆ R + , where θ ∈ Θ ⊆ R q is a set of unknown calibration parameters associated with the compartmental model. For instance, an SIR model contains two parameters, θ = (β, γ), where β is the infection rate controlling the transition from S to I, and γ is the recovery rate controlling the transition from I and R. Suppose that y i is the reported number of infected cases at time x i . Then, given the reported number of infected cases in n days,
, the commonly used approach to estimate the parameters is to minimize the sum of squared differences between actual numbers of infected cases and simulation outputs from compartmental models. The estimated parameters are denoted byθ LS n , where LS stands for least squares, and they are obtained bŷ
(1)
Despite the wide applications of the least squares approach, it can be shown that the least squares estimator does not achieve the semiparametric efficiency when the simulator f (x, θ)
is imperfect, meaning that the simulation output cannot perfectly fit the response, even with the best fit of θ. The asymptotic variance can be reduced by the proposed estimator introduced in this subsection. Theoretical justifications are provided in Section 3.
Assume that the number of cases y i follows a Poisson distribution: y i ∼ Poi(λ(x i )) for i = 1, . . . , n, and y i and y j are mutually independent for any i = j. λ(x i ) is the true mean function of y i . The function λ(x) is often called the true process in the computer experiment literature (Kennedy and O'Hagan, 2001; Wu, 2015, 2016) . Ideally, if the underlying mean function λ(x) is known, the true parameter can be defined as the minimizer of the expected squared discrepancy between the true process and the simulation output, that is,
If we assume that x is uniformly distributed, then it leads to the minimization of the L 2 projection of the discrepancy. That is,
where g L 2 (Ω) = Ω g 2 1/2 .
In reality, the underlying true process λ(·) is unknown that needs to be estimated by observed data. Therefore, given the data {(x i , y i )} n i=1 , we propose to estimate the true process by the kernel Poisson regression (van de Geer, 2000; Shim and Hwang, 2011) . That is,λ n (·) = exp{ξ n (·)} and ξ n = arg min
where · 2 N Φ (Ω) is the norm of the reproducing kernel Hilbert space generated by a given positive definite reproducing kernel Φ, and κ n is a tuning parameter, which can be chosen by cross-validation methods. Thus, the proposed estimator of θ, which we call L 2 -estimator throughout this paper, is the minimizer of the L 2 projection as follows:
The optimal solution of (3) has the form ofξ
can be obtained by the iterative re-weighted least squares algorithm (Green and Yandell, 1985; Hastie and Tibshirani, 1990; Wahba et al., 1995) . The detail of the algorithm is given in Appendix A. In practice, the calculation of the L 2 norm in (4) can be approximated by numerical integration methods, such as Monte Carlo integration (Caflisch, 1998) .
For some stochastic agent-based simulations in epidemiology that can be computationally intensive, it is infeasible to obtain f (x, θ) by conducting simulations for all possible combinations of the input parameters. In such cases, the simulator is often approximated by a computationally cheaper, statistical emulator, which is denoted byf (x, θ), and there are extensive studies on the development of statistical emulators in the computer experiment literature (Santner et al., 2018) . For these simulators, the proposed L 2 -estimator can be similarly obtained by replacing f (x, θ) in (4) with its emulatorf (x, θ). The applications of the proposed method with various existing emulators are demonstrated in Sections 4 and 5.
Theoretical properties of the L 2 -estimator are discussed in this section, including the asymptotic consistency and the semiparametric efficiency. Theoretical comparisons with the least squares estimators are also provided by examining their asymptotic variances. The proofs are given in Supplementary Material.
The following theorem shows that the L 2 -estimatorθ n in (4) is asymptotically consistent and normally distributed.
Theorem 3.1. Under the regularity conditions C1-C10 in Supplementary Material S2, we
By the Delta's method, the following corollary extends the result of Theorem 3.1 to a function of the L 2 -estimator, which we denote as g(θ n ).
Corollary 3.2. For a function g satisfying the property that ∇g(θ * ) exists and is non-zero valued, we have
as n → ∞.
Corollary 3.2 provides a theoretical support for the estimation and inference of some commonly used quantities of interest in epidemiology, such as the basic reproduction rate, which measures the transmission potential of a disease. For instance, in an SIR model the basic reproduction rate is a ratio of two calibration parameters, that is, g(θ 1 , θ 2 ) = θ 1 /θ 2 , where θ 1 = β and θ 2 = γ. The result of Corollary 3.2 can then be applied to construct the confidence intervals for the basic reproduction rate.
When estimating the unknown parameters in compartmental models, the parameters of interest θ are q-dimensional, while the parameter space of model (4) contains an infinite dimensional function space which covers λ. Therefore, the calibration problem is regarded as a semiparametric problem. For these problems, the estimation method that can reach the highest estimation efficiency is called semiparametric efficient (Bickel et al., 1993; . It can be shown in the following theorem that the proposed L 2 -estimator is semiparametric efficient.
Theorem 3.3. Under the regularity conditions in Theorem 3.1,θ n is semiparametric efficient.
Remark 3.4. Theorems 3.1 and 3.3 remain valid even if the simulator f is replaced by its emulatorf , provided that the emulatorf satisfies the conditions C11-C12 given in Supplemental Material S2.
In the next theorem, the asymptotic properties of the least squares estimator are developed and compared with those of the L 2 -estimator.
Theorem 3.5. Under the regularity conditions C1-C4 and C13-C14 in Supplementary
Similar to the L 2 -estimator, it is shown that the least squares estimator is asymptotically consistent and normally distributed. It can also be shown that W 1 ≥ W 0 , which leads to
This implies that the asymptotic variance of the least squares estimatorθ LS n is greater or equal to that ofθ n . The equation holds if and only if
This result indicates that, if ∂f ∂θ (x, θ * ) = 0 for all x ∈ Ω, then (6) holds only if λ(x) = f (x, θ * ) for all x ∈ Ω, which implies that the least squares estimatorθ LS n is less efficient thanθ n if f is an imperfect simulator, i.e., λ(x) = f (x, θ * ) for some x ∈ Ω.
In this section, two artificial examples are conducted to examine the finite sample performance of the proposed method and compare the estimation performance with the least squares approach. These numerical studies are conducted on a desktop with 3.5 GHz CPU and 8GB of RAM, and 4 CPUs are available for parallel computing.
We consider an imperfect simulator adapted from Tuo and Wu (2015) with one calibration parameter. The true process is assumed to be λ(x) = exp(x/2) sin(x/2) + 3, where
x ∈ Ω = [0, 2π], and it is illustrated in the left panel of Figure 1 as the solid line. The data are generated from equal-spaced inputs in [0, 2π] with n = 50 and the outputs are generated from a Poisson distribution with the mean process λ(x i ) for i = 1, . . . , 50, which are shown as the solid dots in the left panel of Figure 1 .
We assume that the simulation output is
This simulator is imperfect because √ θ 2 − θ + 1 is always positive for any θ ∈ Θ. The true parameter can be analytically solved by minimizing (2), which gives that θ * = −0.1789. Plugging in the true calibration parameter, the simulator f (x, θ * ) is demonstrated as the dashed line, which is imperfect because, even with the true minimizer, the discrepancy between the simulation output and the true process is nonzero.
The performance of the L 2 -estimator is compared with the least squares estimator based on the mean squared errors (MSE) obtained from 100 replicates, that is, 100 i=1 (θ i − θ * ) 2 /100, whereθ i is the estimate at the i-th replicate. Their MSEs are shown in the first two bars in the right panel of Figure 1 , which are decomposed into squared biases (dark red) and variances (light blue). It shows that the L 2 -estimator ("L2") provides about 51% smaller MSE than the least squares estimator ("LSE"). Such a reduction is mainly due to the reduction in the estimation variance of the L 2 -estimator. On the other hand, the squared biases are relatively small for both estimators, which is consistent with the asymptotic result of estimation consistency in Section 3.
We further compare the estimation performance for the cases when the simulations are computationally demanding and therefore statistical emulators are built as surrogates. Before comparing the estimation performance, we first examine the emulation performance of two existing emulation methods that are applicable to count data, which are the multiresolution functional ANOVA emulation (Sung et al., 2020b) Table 1 : Emulation performance for the example with one calibration parameter, where m is the sample size of unique locations and a is the number of replicates. RMSPEs are reported for the two emulators based on 10, 000 random predictive locations.
Next, we compare the estimation performance with the hetGP emulator built by m = 100, a = 100 samples. For both estimators, the actual simulator f (x, θ) is replaced by the hetGP emulator. The MSEs of the L 2 -estimator ("L2+emulator") and the least squares estimator ("LSE+emulator") are shown in the last two bars in the right panel of Figure   1 . Similar to the previous result without emulators, the L 2 -estimator provides a smaller MSE, squared bias, and variance compared to the least squares estimator. By comparing the first two and last two bars, it is not surprising to see that the MSEs of "L2+emulator"
and "LSE+emulator" are larger than "L2" and "LSE" due to the prediction error from emulation. It is also worth noting that the second bar is higher than the first and the third, which indicates that, by choosing the proposed estimator, the MSE can be smaller than the least squares estimator even if the simulations are approximated by an emulator.
We consider a more complex problem with three calibration parameters adapted from Plumlee (2017) . Assume that the true mean process is λ(x) = 3x + 3x sin(5x) + 3 and the (2). MSEs. Similar to the previous example, the MSEs are decomposed into squared biases (dark red) and the estimation variance (light blue). From these results, it appears that the L 2 -estimator outperforms the least squares estimator for all of the three parameters by providing a smaller squared bias, variance, which result in a smaller MSE.
In this example, we also examine the prediction performance of the two existing emulators, the unique input locations are sampled from (x, θ) ∈ (χ, Θ) ⊆ R 4 using an LHD. After the emulators are built, the RMSEs are computed based on the predictions of 10, 000 untried input locations, and the prediction results are summarized in Table 2 with different settings of m and a. Similar to the previous example, the hetGP method outperforms MRFA in terms of prediction accuracy and computational time. With a larger a, i.e., more replicates, the prediction accuracy of hetGP appears to increase without much increase in computational time. Thus, we select hetGP as the emulator in the following analysis.
We now compare the estimation performance for the cases where emulators are constructed as surrogates to the actual simulations. The emulator is built by hetGP with m = 300, a = 100 and based on the emulator, the estimation performance is summarized by the last two bars in each of the three plots in Figure 2 . The results indicate that, either when the actual simulator is conducted or emulated, the L 2 -estimator provides smaller squared biases, estimation variances, and MSEs compared to the least squares estimator. 
The proposed method is applied to estimate unknown parameters in compartmental models for the analysis of COVID-19 pandemic. According to the epidemiology literature, the SEIR model, which consists of four compartments, S usceptible-E xposed-I nfectious-Recovered, is widely recommended for COVID-19 simulations because it accounts for the incubation period through the exposed compartment Carcione et al., 2020; Mwalili et al., 2020; He et al., 2020; Annas et al., 2020) . We focus on two types of SEIR simulators: a deterministic simulator and a stochastic simulator, which are discussed in Sections 5.1 and 5.2, respectively. To estimate the unknown parameters in these two SEIR simulators, we collect the actual numbers of infected cases from Johns Hopkins University CCSE repository (Dong et al., 2020) through an R package covid19.analytics (Ponce, 2020) . For each country, there are 107 observations collected from March 1st, 2020, to June 15th, 2020, denoted by y i , where i = 1, · · · , 107. The studies are conducted for the top 20 countries which have the highest cumulative confirmed cases reported on June 15th, 2020.
The SEIR model extends the idea of SIR to further consider an E compartment before entering the I compartment, which indicates a significant incubation period during which individuals have been infected but are not yet infectious themselves. Mathematically, a deterministic SEIR model can be written as:
where S, E, I and R represent the numbers of cases in the corresponding compartment, N = S + E + I + R is the total population, x is time, β is the contact rate that represents the average number of contacts per person per time in the susceptible compartment, γ is the recovery rate from the infectious compartment, and κ is the incubation rate which represents the rate of latent individuals becoming infectious, or equivalently, the average incubation period is 1/κ. These ordinary differential equations have no analytical solution but can be solved by numerical solvers, such as the ODEPACK (Hindmarsh, 1983) . The basic reproduction number, R 0 = β/γ, is the expected number of new infected cases from an infectious individual in a population where all subjects are susceptible. It is of great interest in epidemiology because it provides the insight of the virus transmission, which is essential to the development of effective strategies to contain the pandemic. An accurate and efficient estimation of R 0 is not only important for the public safety, but it also has significant impacts on global economy.
This model includes five unknown parameters: β, κ, γ and the initial numbers of infectious and exposed cases (denoted by I 0 and E 0 , respectively), so we have θ = (β, κ, γ, I 0 , E 0 ). For each country, the L 2 -estimator of θ is obtained by minimizing (4) Plugging in the L 2 -estimators, the simulation results for the top 12 countries with the highest R 0 are demonstrated as the red curves in Figure 4 . In general, it appears that the simulation results can reasonably capture the overall trend observed from the actual numbers of infected cases, which are shown as the gray dots. For Canada, France, and Turkey, the discrepancy between the simulation results and actual observations is relatively larger than the other countries. This is partly because SEIR is an imperfect simulator which is built based on some assumptions or simplifications, and these assumptions may have larger deviations from the reality for certain countries. Another reason is that the intrinsic dynamics are neglected in the deterministic simulations. To take into account the dynamics, a stochastic simulator is considered in the next subsection.
A stochastic SEIR simulation provides a more sophisticated and realistic framework to integrate infection dynamics in different compartments as continuous-time Markov chains (Allen, 2008; Andersson and Britton, 2012; Allen, 2017) . To conduct these simulations, we implement an R package, SimInf (Widgren et al., 2019) , in which the simulation results are obtained by the Gillespie stochastic algorithm (Gillespie, 1977) . Stochastic SEIR simulations are computationally more demanding. For example, it takes more than 10 minutes to produce one simulation result for one country under a given parameter setting.
It is computationally infeasible to perform simulations for all the possible combinations of the parameters; therefore, an emulator is constructed as an efficient surrogate to the actual simulation. In this study, we consider the hetGP emulator, which is built based on the simulations generated using a 60-run LHD for parameter settings with 19 equal-spaced time steps in x, which leads to the total sample size of m = 1140. For each parameter-input setting, 100 replicates are simulated, i.e., a = 100. Based on this emulator, it takes less than two seconds to emulate the result for an untried parameter setting, which is significantly faster than the actual stochastic simulation.
By replacing the actual simulator with the hetGP emulator, the L 2 -estimators are obtained by minimizing (4). The corresponding estimates of R 0 and their 95% confidence intervals are summarized in Figure China and Germany remain the smallest which agrees with the findings in the deterministic simulations, but Germany appears to have higher R 0 than the one obtained from the deterministic simulations. We further report the estimated incubation period, 1/κ, for each country and the corresponding 95% confidence intervals in Figure 6 . The overall average incubation period is 6.31 as indicated by the red dashed line.
In Figure 7 , the actual numbers of infected cases are illustrated as the gray dots. By plugging in the L 2 -estimators, the simulation results for the top-12 countries with the highest R 0 are illustrated as the red curves, along with the 95% predictive intervals as the red-shaded regions. Overall, the simulation results show a much better agreement with the actual observations compared to the deterministic ones in Section 5.1. In particular, by taking into account the intrinsic dynamics, the simulation discrepancy for Canada is significantly reduced from the deterministic one shown in Figure 4 . Apart from the frequentist approach studied in this paper, we are currently developing a Bayesian framework that extends the recent developments of Bayesian calibration to count data. For example, the orthogonal Gaussian process models (Plumlee et al., 2016) or the projected Gaussian process (Tuo, 2019) can be used to model the model discrepancy, which addresses the unidentifiability issue for continuous outputs, and it is conceivable to further extend the modeling to count data by incorporating the idea of the generalized calibration in Grosskopf et al. (2020) . Another interesting direction that deserves further studies is to relax the constant parameter assumption. Instead, the calibration parameters can be assumed to be functions of some factors, such as time or temperature, which not only increases the model flexibility but also can provide further insights to the time-course dynamics of the COVID-19 infection.
Acknowledgements: This work was supported by NSF DMS 1660477 and NSF HDR TRIPODS award CCF 1934924.
A Algorithm to Estimateξ n in (3) Since the optimal solution has the form of ξ n (x) = b + n i=1 a i Φ(x i , x), one can show that the penalized likelihood in (3) can be rewritten as
where a = (a 1 , . . . , a n ), φ(x) = (Φ(x, x 1 ), . . . , Φ(x, x n )), and Φ = (Φ(x i , x j )) 1≤i,j≤n . The optimal solution of a and b can then be obtained by taking the first-order partial derivatives of the objective function with respect to a and b and setting them equal to zero, which can be solved by the iterative re-weighted least squares algorithm as follows. Denote
where 1 n = [1, . . . , 1] T and 0 n = [0, . . . , 0] T , and denote W as an n × n diagonal matrix with diagonal elements W ii = exp(b + a T φ(x i )). Then, in each step, one first solve for
with an initial guess of η, which is a vector of size n, and then update each element of η by
The estimateβ can then be obtained by continuing solving for β and η iteratively until some convergence criterion is met. Supplementary Materials for "Efficient Calibration for Imperfect Epidemic Models with Applications to the Analysis of COVID-19" S1 Asymptotic Results ofξ n of (3)
We start with a result developed by van de Geer (2000) for general nonparametric regression (Lemma 11.4 and 11.5 in van de Geer (2000)). Denote
Suppose F is the class of all regression functions equipped with the Sobolev norm · H m (Ω) , which is defined by
. Letξ n := arg min
for some κ n > 0. Then the convergence rate ofξ n is given in the following lemma.
Lemma S1.1. Let ξ 0 ∈ F. Assume that there exists some nonnegative k 0 and k 1 so that
where z 0 = ξ 0 (x) and x ∈ Ω. For κ −1 n = O(n 2m/(2m+1) ), we have
In fact, the norms of some reproducing kernel Hilbert spaces (RKHS) are equivalent to Sobolev norms. For instance, the RKHS generated by the Matérn kernel function, given by
where ν ≥ 1 and ρ ∈ R + are tuning parameters and K ν is a Bessel function with parameter ν, is equal to the (fractional) Sobolev space H ν+d/2 (Ω), and the corresponding norms · N Φ (Ω)
and · H ν+d/2 (Ω) are equivalent (Wendland, 2004; Tuo and Wu, 2016) . Therefore, as a consequence of Lemma S1.1, we have the following proposition forξ n obtained by (3).
Proposition S1.2. Suppose that ξ 0 ∈ F = N Φ (Ω), and N Φ (Ω) can be embedded into H m (Ω). Then, for κ −1 n = O(n 2m/(2m+d) ), the estimatorξ n in (3) andλ n = exp{ξ n } satisfy
and ξ n − ξ 0 L 2 (Ω) = O p (κ 1/2 n ).
Proposition S1.2 suggests that one may choose κ n n −2m/(2m+d) to obtain the best convergence rate λ n − λ 0 L 2 (Ω) = O p (n −m/(2m+d) ), where a n b n denotes that the two positive sequences a n and b n have the same order of magnitude.
S2 Regularity conditions for Theorems 3.1, 3.3, 3.5, and Remark 3.4
The regularity conditions for Theorem 3.1 are given below. For any θ ∈ Θ ⊂ R q , write θ = (θ 1 , . . . , θ q ). Denote e i = y i − λ(x i ) and ξ(·) = log λ(·). C2: θ * is the unique solution to (2) and is an interior point of Θ.
C3: sup θ∈Θ f (·, θ) L 2 (Ω) < +∞.
C5: There exists a neighborhood U ⊂ Θ of θ * such that
for all θ ∈ U and all i, j = 1, . . . , q.
The Donsker property is an important concept in the theoretical studies of empirical processes. The definition and detailed discussion are referred to van der Vaart and Wellner (1996) and . Tuo and Wu (2015) showed that if the conditions of Proposition S1.2 hold and m > d/2, then N Φ (Ω, ρ) is a Donsker. The authors also mentioned that under the condition C1 and E[exp{C|y i − λ(x i )|}] < +∞ for some C > 0, the conditions of Proposition S1.2 are satisfied. Therefore, by choosing a suitable sequence of κ n , say κ n n −2m/(2m+d) , one can show that condition C10 holds and C8 and C9 are ensured by Proposition S1.2.
Additional regularity conditions for Theorem 3.4 are given below, where we assume that there are more than n samples in the computer experiments for building the emulatorf .
Additional regularity conditions for Theorem 3.5 are given below.
C13: Suppose that there exists a neighborhood U of θ * such that f (x, ·) ∈ C 2,1 (U ) for all x ∈ Ω, where C 2,1 denotes the space of functions whose second derivatives are Lipschitz.
C14:
In this section, we replace the simulator f in (4) We first show that under the regularity conditions C1-C12 in Supplementary Material S2,θ
The proof is developed along the lines described in Theorem 1 of Tuo and Wu (2015) and
Theorem 3.1 of Sung et al. (2020) . We first prove the consistency,θ n p − → θ * . It suffices to prove that λ n (·) −f (·, θ) L 2 (Ω) converges to λ(·) − f (·, θ) L 2 (Ω) uniformly with respect to θ ∈ Θ in probability, which is ensured by
where the inequality follow from the Schwarz inequality and the triangle inequality.
Since it can be shown that f L 2 (Ω) ≤ Vol(Ω) f L∞(Ω) , for any f ∈ L ∞ (Ω), where Vol(Ω) is the volume of Ω, we have
By combining (S3.5) and (S3.6) and the conditions C8, C9, and C11, we have that (S3.4) converges to 0 uniformly with respect to θ ∈ Θ, which proves the consistency ofθ n .
Sinceθ n minimizes (S3.3), by invoking conditions C1, C2 and C5, we have
and by the conditions C8, C11 and C12, it implies (3), we know thatξ n minimizes l over N Φ (Ω). Sinceθ n p − → θ * and by condition C5, ∂f ∂θ j (·,θ n ) ∈ N Φ (Ω) for j = 1, . . . , q with sufficiently large n. Then we have
is Donsker. Thus, by Theorem 2.10.6 in van der Vaart and Wellner (1996) , F 1 = {exp{g} − exp{ξ} : g ∈ N Φ (Ω, ρ)} is also Donsker because the exponential function is a Lipschitz function. By condition C5, the class F 2 = { ∂f ∂θ j (·,θ n ), θ ∈ U } is Donsker. Since both F 1 and F 2 are uniformly bounded, by Example 2.10.8 in van der Vaart and Wellner (1996) the product class F 1 × F 2 is also Donsker. Thus, the asymptotic equicontinuity property holds, which implies that for any > 0 there exists a δ > 0 such that lim sup n→∞ Pr sup Suppose ε > 0 is a fixed value. Condition C9 implies that there exists ρ 0 > 0 such that P r( ξ N Φ > ρ 0 ) ≤ ε/3. In addition, choose δ 0 to be a possible value of δ which satisfies (S3.9) with = ε/3 and ρ = ρ 0 . Condition C8 implies that Pr( exp{ξ n }−exp{ξ} L 2 (Ω) > δ 0 ) < ε/3. Then, for sufficiently large n, we have
The first and third inequalities follow from the definition ofξ • n , and the last inequality follows from (S3.9). Thus, this implies that E 1n (ξ n , θ) tends to zero in probability, which
Then, by substituting (S3.7) to (S3.10) and using condition C2, Taylor expansion can be applied to (S3.10) at θ * , which leads to
whereθ n lies betweenθ n and θ * . By the consistency ofθ n , we then haveθ n p − → θ * , which implies that
Thus, we have
Next, we consider D n . Define the empirical process
where e i = y i − λ(x i ) θ ∈ U . Condition C1 implies that the set {ζ θ ∈ C(R × Ω) : ζ θ (e, x) = e ∂f ∂θ j (x, θ) − e ∂f ∂θ j (x, θ * ), θ ∈ U } is a Donsker class, which ensures that E 2n (·) converges weakly in L ∞ (U ) to a tight Gaussian process, denoted by G(·). Without loss of generality, we assume G(·) has continuous sample paths. Then, by the continuous mapping theorem (van der Vaart, 1998) and the consistency ofθ n , we have E 2n (θ n ) p − → G(θ * ). Because E 2n (θ * ) = 0 for all n, we have G(θ * ) = 0. Then, we have E 2n (θ) p − → 0, which gives
Lastly, we consider E n . Applying conditions C5, C9, C10, we have
By combining (S3.8), (S3.11), (S3.12) and (S3.13), we havê
Thus, given condition C6 and
(S3.14) implies the asymptotic result of Theorem 3.1 by the central limit theorem, that is,
If suffices to show thatθ n given in (4) has the same asymptotic variance as the estimator obtained by using maximum likelihood (ML) method. Consider the following q-dimensional parametric model indexed by γ, ξ γ (·) = ξ(·) + γ T ∂f ∂θ (·, θ * ), (S4.15) with γ ∈ R q . By the fact that ξ(·) = log λ(·) and (S4.15), it becomes a Poisson regression model with the coefficient γ. The log-likelihood of γ is
and its first and second derivatives are ∂l(γ) ∂γ T = 1 n n i=1 y i ∂f ∂θ (x i , θ * ) − exp ξ(·) + γ T ∂f ∂θ (·, θ * ) ∂f ∂θ (·, θ * ) and ∂ 2 l(γ) ∂γ∂γ T = − 1 n n i=1 exp ξ(·) + γ T ∂f ∂θ (·, θ * ) ∂f ∂θ (·, θ * ) ∂f ∂θ T (·, θ * ) .
Since the true model is y i ∼ Poi(λ(x i )), the true value of γ is 0. Hence, under the regularity conditions of Theorem 3.1, the ML estimator of γ has the asymptotic expression
where W 0 is defined in (5). Then, a natural estimator for θ * in (2) for all t near zero, and let
(S4.17) implies that φ(θ(t), t) = 0 for all t near zero. Then, by the implicit function theorem, we have
By the delta method, we havê θ n − θ * = θ(γ n ) − θ(0) = ∂θ(t) ∂t T t=0γ n + o p (n −1/2 ), and together with (S4.16) and (S4.18), it yieldŝ
Thus, the asymptotic variance of the ML estimator is 4V −1 W 0 V −1 , which is the same as that ofθ n (see Theorem 3.1). Therefore, the L 2 estimatorθ n in (4) is semiparametric efficient.
Consider a parametric model, y i ind.
∼ N (f (x i , θ), 1), which is misspecified because the true model is y i ind.
∼ Poi(λ(x i )). It is easy to see thatθ LS n is the maximum likelihood estimator (MLE) of θ under this misspecified model. Thus, under the regularity conditions C1-C4 and C11, the MLE converges to θ in probability which uniquely minimizes Kullback-Liebler divergence (Huber, 1967; . That is,
where h 0 and h 1 are the density functions of y i under the misspecified model and the true model, respectively, and E 0 denotes the expectation with respect to y i ∼ h 0 . Since the true model follows a Poisson distribution, we have
By the law of large numbers, we have
Therefore,θ LS n p → θ = arg min θ∈Θ E[λ(X) − f (X, θ)] 2 , which shows thatθ LS n is consistent because θ = θ * .
In addition, by , the asymptotic normality ofθ LS n can be shown as follows, (S5.20) and 