The Interplay between COVID-19 and the Economy in Canada

: We propose a generalized susceptible-exposed-infected-removed (SEIR) model to track COVID-19 in Canadian provinces, taking into account the impact of the pandemics on unemployment. The model is based on a network representing provinces, where the contact between individuals from different locations is deﬁned by a data-driven mixing matrix. Moreover, we use time-dependent parameters to account for the dynamical evolution of the disease incidence, as well as changes in the rates of hospitalization, intensive care unit (ICU) admission, and death. Unemployment is accounted for as a reduction in the social interaction, which translates into smaller transmission parameters. Conversely, the model assumes that higher proportions of infected individuals reduce overall economic activity and therefore increase unemployment. We tested the model using publicly available sources and found that it is able to reproduce the reported data with remarkable in-sample accuracy. We also tested the model’s ability to make short-term out-of-sample forecasts and found it very satisfactory, except in periods of rapid changes in behavior. Finally, we present long-term predictions for both epidemiological and economic variables under several future vaccination scenarios.


Introduction
Since the seminal work of Kermack and McKendrick (1927), mathematical models have played an essential role in the modeling of infectious diseases dynamics. Susceptibleinfected-removed (SIR) and susceptible-exposed-infected-removed (SEIR) models are two of the most widely used models designed to describe the spread of an emerging pathogen in a susceptible population (Keeling and Rohani 2008), and since the beginning of the coronavirus disease 2019 (COVID-19) pandemic, several SEIR-type models have been proposed to describe the spread of the SARS-CoV-2 virus Calvetti et al. 2020b;Campos et al. 2021;Chen et al. 2022;Gatto et al. 2020;Grimm et al. 2021;Lyra et al. 2020).
During the first outbreak of COVID-19, several catastrophic predictions (regarding the number of deaths and hospital occupations) have emerged worldwide (Ioannidis et al. 2020). Such predictions were based, in general, on over-simplistic models and poor data. Fortunately, such forecasted scenarios did not materialize, although more accurate predictions were urgently needed. To address the need of better data, public authorities from different countries have improved their data collection, data organization, and data sharing services on, among others, COVID-19 testing, infections, deaths, and hospitalizations. This allowed the World Health Organization (WHO) to track the worldwide evolution of the COVID-19 pandemics and let researchers test their models with real data 1 .
To improve the forecast ability of SEIR-type models, a series of features were incorporated into the original models. One of them, inspired by the widely implemented mobility restrictions, was the spatial spread of the pathogen. Some works considered the spatial dynamics based on the distance between neighborhoods Calvetti et al. 2020b;Gatto et al. 2020;Wu et al. 2020), others used mobility data based on mobile phone information (Alleman et al. 2021;Candido et al. 2020;Coelho et al. 2020;Kishore et al. 2021). Another important feature incorporated into the model was the different levels of disease severity, commonly grouped as asymptomatic, mild, severe, and critical Campos et al. 2021). Such severity was highly dependent on demographic factors, such as age and sex (Bhopal and Bhopal 2020;CDC 2020;Jin et al. 2020;Verity et al. 2020;WHO 2020), which were also considered in some models Cartocci et al. 2021;Grimm et al. 2021), given their importance in the prediction of hospital and intensive care unit (ICU) usage. The performance and forecast accuracy of such SEIR-type models were also considerably improved when some parameters were considered as time-dependent (Albani et al. , 2022Calvetti et al. 2020aCalvetti et al. , 2020b. This allowed the models to incorporate regime changes in the dynamics caused, for example, by the emergence of new variants or treatments that change the pathogen transmission pattern and mortality rates, as well as a variety of public health measures, ranging from lockdowns to vaccination. These modeling improvements resulted in a series of theoretical and data-driven numerical experiments to quantitatively access the performance of different techniques used to control COVID-19 outbreaks. For example, SEIR-type models were used to quantify the efficacy of different non-pharmaceutical measures, such as lockdowns and the use of masks (Alleman et al. 2021;Eikenberry et al. 2020; IHME COVID-19 Forecasting Team 2020; Kishore et al. 2021). Epidemiological models were also used to predict the impact of vaccination strategies in the population (Albani et al. 2021b;Libotte et al. 2020), as well as to evaluate the effects of underreporting of infections and deaths (Albani et al. 2021a;Aronna et al. 2022;Chen et al. 2022;Kupek 2021;McMahan et al. 2021;Saberi et al. 2020). It was also possible to access the role of socioeconomic programs in the mitigation of outbreaks .
Concerning the calibration of epidemiological models, different techniques were applied, including deterministic (optimization) and stochastic (Bayesian inference) methods. The main issue in calibration is overfitting, where a model fits observed data well, but has serious difficulty in predicting future scenarios. This phenomenon is caused by the interplay between uncertainties and noise associated with the data collection and the use of an excessive number of model parameters. To address this problem, different techniques can be applied depending on the calibration method. In general, deterministic models use some sort of regularization, for example, à la Tikhonov (Engl et al. 1996;Scherzer et al. 2008). In the context of stochastic methods, inference techniques based on statistical analysis are used to incorporate such uncertainties into the solution (Gamerman and Lopes 2006;Migon et al. 2014;Somersalo and Kapio 2004), at the expense of being more computationally intensive than deterministic methods. In the context of the calibration of epidemiological models for COVID-19, deterministic methods were used, for example, by ; Chen et al. (2022); He et al. (2020); Leonov et al. (2021); Ramadijanti and Basuki (2020), and stochastic techniques were used in Calvetti et al. (2020aCalvetti et al. ( , 2020b Gianfelice et al. (2022); Lyra et al. (2020); Monod et al. (2021); Unwin et al. (2020).
The economic impacts of the COVID-19 pandemic were also widely studied under different perspectives. Some studies (Fernández-Villaverde and Jones 2020; Maliszewska et al. 2020;OECD 2020) compiled and summarized the evidence of the impact on traditional economic measures, such as GDP, employment, and trade, whereas others (Ashraf 2020;Storm 2021) used statistical methods to quantify the effects of different non-pharmaceutical interventions and compare their effectiveness across multiple countries. As informative as these empirically-based studies can be, one difficulty with this approach is the lack of an underlying model governing the dynamics of both the pandemic and the economy, which makes it particularly challenging for them to make reliable medium-and longterm forecasts.
A different strand of the literature, represented for example by Acemoglu et al. (2021); Alvarez et al. (2021); Bayraktar et al. (2021) focused instead on the design of optimal policies, typically by extending an SEIR-like model through the introduction of an optimal control problem. This type of approach, however, faces the difficulty of unrealistic modeling assumptions, such as the definition of an agreed-upon objective function for a central planner, with well-understood effects of the pandemic on the economy, and the ability to execute precise controls, such as having exact time-dependent fractions of the population under a lockdown regime. A related strand, represented for example by Eichenbaum et al. (2021);Farboodi et al. (2021), seeks to endogenize behavior, such as physical distancing by having rational individuals decide on their level of social interaction based on the probability of becoming infected. As in traditional macroeconomic models, these decisions are made with the view of maximizing the utility of consumption over an infinite time horizon, with the resulting model suffering from the same deficiencies as mainstream dynamic stochastic general equilibrium (DSGE) models, including a lack of heterogeneity, perfect foresight, and instability of the solution path.
As an alternative to either a purely empirical approach or an unrealistic theoretical approach, the work presented, for example, in Kemp-Benedict (2020), Cruz-Aponte and Caraballo-Cueto (2021), and Balmford et al. (2020), combined an epidemiological SEIR-like model with a more phenomenological macroeconomic model. Along similar lines, we highlight the contribution of Góes and Gallo (2021), to which our paper is closely related. Their central idea is to relate a key economic variable, namely the unemployment rate u, with one key epidemiological variable, namely the positivity rate i, through a stylized model of the form The motivation for this model is that, absent the pandemic, the unemployment rate converges to a long-term equilibrium u n , whereas an increase in the number of infected people leads to higher unemployment through the interaction term ρu(t)i(t) in Equation (1) with a constant ρ > 0. This is meant to include all the effects of the pandemic on the economy, both direct and indirect, such as sick people not being able to work, slowdowns in economic activity caused by lockdown measures, business closures, and the like. Conversely, an increase in unemployment rate leads to a decrease in the number of infected people through the term −γu(t) in (2) with a constant γ > 0. Again, this is meant to encapsulate all the effects of the economy on the pandemic, both direct and indirect, such as lower mobility when unemployed people stay at home, as well as overall lower contact and transmission during lockdown periods and corresponding lower economic activity. The equilibrium points for (1) and (2) are namely, an endemic equilibrium with u 1 > u n and a non-endemic one with u 2 = u n . Góes and Gallo (2021) then chose the parameter β based on estimates available in the literature for the contact, recovery, and death rates for COVID-19 in the US, and estimated the parameter γ based on the observed reduction in the basic reproductive number when containment measures were introduced. The parameters α and ρ were estimated using the time series for the unemployment and positivity rates for the US, whereas the parameter u n was chosen to be compatible with conventional estimates of long-term unemployment for the US economy. They simulated the model with these parameter values and observed convergence to the endemic equilibrium with unemployment close to 10% and a positivity rate close to 4.4%. Finally, Góes and Gallo (2021) simulated a stochastic version of (1) and (2) and compared it with the observed dynamics for these two variables, obtaining relatively good agreement from the start of the pandemic until November 2020, using data available prior to the acceptance date of the article in December 2020.
Our work can be characterized as an extension of the model of Góes and Gallo (2021) applied to the Canadian economy. On the epidemiological side, we replaced Equation (1) with a fully developed SEIR-type model with each compartment stratified by location, age, and employment status. On the economic side, we replaced Equation (2) with equations describing the flows between the employed and unemployed groups in each compartment. In addition, similar to Góes and Gallo (2021), we estimated a relationship between the unemployment rate and economic output, known as Okun's law, before and throughout the pandemic, in order to convert the labor market information into predictions of GDP for each province in Canada.
More specifically, in comparison to the existing literature, the main contributions of the present article are the following: 1.
An epidemiological model that accounts for spatial spread, time-dependent parameters, different levels of disease severity, vaccination, age-range dependency, and labor status of those of working age.

2.
Model parameters that are calibrated from the daily Canadian provincial reports on COVID-19 infections, hospitalizations, admissions to ICU, deaths, and monthly labor status.

3.
An explicit connection between the pandemic and the economy through labor market dynamics.
The calibrated model is able to provide accurate predictions and presents a remarkable fit to the observed data, including periods of rapidly varying policies, both economic and epidemiological. Toward the end of the article, we present long-term predictions for the future evolution of COVID-19 and its economic impact based on a number of scenarios for vaccination and loss of immunity.

The Model
We present the model under a general formulation considering different locations , age groups a and employment status s. For the empirical application, = 1, . . . 13 will correspond to different Canadian provinces and territories, a = 1, 2, 3 to the age groups 0 to 19, 20 to 69, and over 69 years old, and s = e, u corresponding to employed and unemployed individuals.
Thus, for each location and age range, the population is distributed into the following compartments: susceptible (S s ), vaccinated (V s ), exposed but not infective (E s ), infective and mildly ill (I s M ), infected and hospitalized (I s H ), infected and in an intensive care unit (I s I ), and recovered (R s ). The deceased individuals are denoted by D, regardless if they came from the employed or unemployed group. Before describing the disease dynamics in more detail, let us introduce a vector notation as follows. For each employment status s, let S s,a, denote the number of susceptible individuals in the age range a and location and define where A denotes the number of age groups and L denotes the number of different locations. That is to say, S s is an N × 1-row vector with N = A × L components, whose first A components are the number of susceptible individuals of different ages in the first location, the next A components are the number of susceptible individuals of different ages in the second location, and so forth. The vectors V s , E s , I s M , I s H , I s I , R s , and D s are defined similarly. We also introduce the tensor product The schematic representation of the model for each location and age range in Equations (7)-(21) can be found in Figure 1.
As in a regular SEIR model, susceptible individuals can become exposed upon contact with infective ones, whereas exposed individuals become infective after an incubation time. Each mildly ill infective individual can recover, become hospitalized, or die. Hospitalized infected individuals can recover or develop more severe symptoms and be admitted to an ICU. Those individuals in the ICU can only recover or die. We also admit the possibility of reinfection, i.e., individuals in the vaccinated and recovered compartments can become susceptible after a while. In principle, individuals in each compartment could move between employment statuses, but we assume that there is no possibility of becoming employed or unemployed while in the hospital. That is to say, only the individuals in the compartments (i.e., susceptible, vaccinated, exposed, mildly infective, and recovered) can change the employment status.
Accordingly, the detailed epidemiological model is defined as follows: The matrices β M , β H , and β I are the transmission parameters for the mildly ill, hospitalized, and in intensive care infective classes, respectively, taking into account the age group and location for each class. More specifically, β M , β H , and β I are time-dependent N × N symmetric matrices, where N = A × L, A is the number of age groups and L is the number of locations. The matrix β M is defined as, where C i,j , for i, j = 1, . . . , L, are A × A symmetric matrices and β i (t), for i = 1, . . . , L, are time-dependent scalar quantities. Notice that the entries for the A × A matrices C i,j could, in principle, be estimated from information about transmission between different age groups in different provinces. Because this is not available in the data, we make the simplifying assumption these entries are identical for a given C i,j , which is equivalent to assuming that, given a pair of provinces (i, j), the probability of transmission is uniform across age groups. Specifically, for each i = 1, . . . , L, all the entries of C i,j are given by with b i,i = 1, and b i,j is the probability of an individual from place i to meet someone from place j. The parameters b i,j depend on the distance between the locations i and j, and are estimated from the daily reports of infections. Again, to reduce the number of parameters, we assume that b i,j can only take five different values b 1 , . . . , b 5 when i = j, depending on whether the distance between the centers 2 of i and j falls is less than 1000 Km, from 1000 to 2000 Km, from 2000 to 3000 Km, from 3000 to 4000 Km, or larger than 4000 Km 3 . The matrices β H and β I satisfy β H = aβ M and β I = bβ M , where a and b are scale parameters that reduce the transmission rate amongst the hospitalized and in ICU individuals.
Most of the remaining parameters in the model are N × 1 column vectors corresponding to specific transition rates between compartments for each age group and location. Namely, the vectors γ M and γ H are the hospitalization and ICU admission rates, respectively. The vectors ν M , ν H , and ν I are the recovery rates of mildly ill, hospitalized, and intensive care-infected individuals, respectively. The vectors µ M and µ I are the mortality rates of those outside a hospital and in the ICU, respectively. Notice that we assume that individuals outside a hospital or in a critical state (in an ICU) can die since in Canada a number of deaths occurred outside hospitals as pointed in WHO (2020). The vector ν denotes the vaccination rate. The vector σ has components equal to the inverse of the mean time of incubation for each age group and location. Similarly, the vector r has components equal to the inverse of the mean time that someone vaccinated or recovered from the disease becomes susceptible again, for each age group and location. Finally, the vectors λ e and λ u denote, respectively, the transition rates from unemployed to employed status and vice-versa. Individuals outside the working age receive the label e and, in these cases, the corresponding components for the rates λ e and λ u are set to 0.
We also assume that unemployed individuals have a smaller risk of getting infected, and this is represented by the scalar parameter 0 < ε < 1, with ε = 0 corresponding to the case of unemployed individuals having no risk of becoming infected. This means that, when someone is unemployed, it is expected that they will take public transportation less often, as they are not going to work, and will spend more time at home, which means a sort of self-isolation. Since hospitalized individuals are isolated, there is no difference in interacting with employed or unemployed individuals in the hospital or ICU.
Regarding the labor market, observe that the total number of employed and unemployed individuals in the population can be obtained from the other variables of the model as follows: Similarly to Góes and Gallo (2021), we assume that these variables obey the following dynamics: In other words, we assume that the vectors λ e and λ u in the epidemiological model in Equations (7)-(21), are given by λ e =λ e + P S ρ e : (σE e + σE u ) and λ u =λ u + P S ρ u : (σE e + σE u ), where the parametersλ e andλ u represent the baseline transition rates into the employed and unemployed groups, whereas the parameters ρ e and ρ u represent the additional sensitivity of these rates with respect to new infections. The parameter P S is the Canadian population size in 2020, i.e., 38.0 million. This parameter is used to scale up the new infections represented by σE s , s = u, e, in comparison to the proportions of employed and unemployed individuals that have higher magnitudes.

Estimation Procedure
We begin with the model parameters in Equations (7)-(21) that can be obtained directly from available data without resourcing to any sophisticated estimation procedure, namely Starting with the vector parameters, the rates of hospitalization γ M , ICU admission γ H , and mortality rate µ M are treated as time-dependent and estimated according to the available information for each Canadian province, whereas the remaining parameters in (29) are treated as constant and calibrated to values from the related literature.
More precisely, γ M (t) is assumed to be identical for all age groups and labor statuses. It is set such that γ M (t) ∑ a,u I ,a,u M (t) is equal to the number of hospitalizations on day t. The rate γ H (t) is defined similarly, using the daily reports of ICU admissions and hospitalizations in the -th place. The death rate µ I is set to 0.39 for all locations and age groups, following Abate et al. (2020). The death rate µ M (t) is also assumed to be identical for all age groups and labor status. It is defined to be such that µ M (t) ∑ a,u I ,a,u is equal to the reported deaths on day t, in the -th place, minus µ M (t) ∑ a,u I ,a,u M (t), for all age groups. The rates of recovery of mildly ill, hospitalized, and intensive care individuals are set, respectively, as ν M = 1/14 days −1 (WHO 2020), ν H = 1/12 days −1 (Guan et al. 2020), and ν I = 1/9 days −1 (Grasselli et al. 2020) for all locations and age groups. The incubation time 1/σ is set to be constant and equal to 5.1 days for all age groups and locations, which is the median time found by Lauer et al. (2020). Similarly, the rate r is also assumed to be constant and equal for all age groups and locations. For the baseline calibration of the model, we set r = 0, whereas we use chose different values for the scenario analysis presented later. It is now well documented that vaccinated and recovered individuals can lose their immunity, and seroprevalence decreases with time. However, this time can vary dramatically depending on a series of factors (Goldberg et al. 2021;Levin et al. 2021;Spellberg et al. 2021). In the simulations that follow, setting r = 1/360 is equivalent to assuming that 360 days is the mean time vaccinated and recovered individuals return to the susceptible condition.
Finally, the vector ν denotes the vaccination rate and is proportional to the number of individuals receiving the second dose daily in each age group and location times 0.7, which represents the efficacy of the vaccines used in Canada. Other ways of including spatial information into SEIR-like models can be found in Calvetti et al. 2020b;Gatto et al. 2020;Keeling and Rohani 2008).
Moving on to scalar parameters, similarly to , the quantities a and b are set to a = 0.10 and b = 0.01. As we do not have epidemiological information accounting unemployment, for simplicity, we set ε = 0.5, which means that transmission among unemployed individuals is half the transmission among employed people, due to the self-isolation caused by not going to work. Whereas small perturbations around this value for ε do not affect the qualitative behavior of the model in a significant way, a systematic sensitivity analysis of the model for different values of this parameter is not feasible, as it requires the repeated re-estimation and simulation of a high-dimensional model with large computational costs. In the event that disease incidence data for employed and unemployed groups become available, the parameter ε could be estimated similarly to the other parameters in the model.
The remaining model parameters need to be estimated from publicly available daily numbers of COVID-19 infections through a so-called maximum a posteriori technique from Bayesian Inference Somersalo and Kapio (2004). We denote by Θ the following vector of unknown parameters: where I M (0) is the initial population of mildly infective individuals in the th location, regardless of age group and employment status, and c and b 1 , . . . , b 5 are the parameters introduced in (23). We then use available data on newly infected individuals to estimate the vector Θ as follows. Let us assume that the daily number of new infections in the -th location, denoted by I , is Poisson-distributed with parameter σE (t), where E (t) denotes the number of exposed and not infective individuals in the -th location, regardless of age group or employment status. More precisely, using the notation introduced in (5), we have that The logarithm of the likelihood function 4 is then with T the number of elements in the time series. We use non-overlapping consecutive intervals of 20 daily observations each to estimate the parameters in a time-dependent manner.
Regarding the initial conditions for the model, notice that, using the notation introduced in (5) Once the I M (0) parameters are estimated as explained above, we obtain each of the I s,a, M (0) based on the proportion of the population for each age range and labor status. The remaining initial conditions for the model are set as follows: S s,a, (0) is set to the proportion of the Canadian population in the th location at the ath age minus the sum of the initial conditions for the other compartments with the same values for s,a, and . If a = 1 or 3, S u,a, (0) = 0. If a = 2, then S s,2, (0) is multiplied by the proportions of unemployed and employed individuals. Apart from the infective compartments, the initial conditions for all the other compartments are set to zero.
The estimation of the parameters λ e and λ u are performed similarly, but using data about changes in employment status instead. Firstly, we have to adjust the labor market data to be used within the epidemiological model since for each location we have monthly totals of employed and unemployed individuals and daily reports of new infections. As the model does not consider natural births and deaths, which means that the sum of all the compartments is constant in time, we divide the monthly numbers of employed and unemployed by the corresponding monthly labor force (i.e., the sum of employed and unemployed individuals). The resulting numbers are then multiplied by the proportion of individuals of working age in the same location for the year 2020. Moreover, to transform the monthly reports into daily data, we interpolate them linearly. The resulting series are then used to estimate the dailyλ s (t) and ρ s (t), for s = u, e and = 1, . . . , L. Since the number of unknowns for each day is four and the daily data are two, a moving window of four consecutive days is used. For such four consecutive days, we assume that the parameter values are constant in time. If Em and Un represent the adjusted daily reports and Em and Un denote the solution of the model in Equations (26) and (27), using the daily reported infections for the th place, and using as initial condition the adjusted reported employed and unemployed numbers, we obtain the parameters (λ e ,λ u , ρ e , ρ u )(t j+1 ) by minimizing the functional for j = 4, . . . , T − 1 and = 1, . . . , L, where the parameters (λ e ,λ u , ρ e , ρ u )(t j ) are estimated in the previous step 5 .

Numerical Results
In this section, we present the results for the calibration of the model in Equations (7)-(21) using real data. We also perform backtesting by considering the accuracy of short-term model predictions.

Model Calibration
We use the 7-day moving average version of the time series of daily reported infections, hospitalizations, ICU admissions, deaths, and second vaccine doses inoculated from the 13 Canadian provinces and territories, provided by Esri Canada. 6 We also use demographic data provided by Statistics Canada. 7 The labor data consists of the provincial monthly total number of unemployed individuals also provided by Statistics Canada. 8 The time series considered in this numerical example correspond to the period from 19 February 2020 to 28 February 2022.
Starting with the epidemiological parameters, Figure 2 shows the values for the 7-day moving average of the time-dependent parameters for the 6 largest Canadian provinces. As expected, we can observe increases in the values of all the parameters during the successive waves of infections or large outbreaks, as well as other patterns broadly consistent with reported news about the pandemic. For example, while the hospitalization rate exhibited a large peak during the omicron wave in the beginning of 2022, the peaks in death and ICU admissions rates for the same period were similar to those during previous waves of infections. This indicates that omicron caused less severe outcomes than the earlier variants. 9 Moving on to the labor market parameters, Figure 3 presents estimated values of the time-dependent parametersλ e ,λ u , ρ e , and ρ u for the six largest Canadian provinces. We can observe in the figure that all six provinces experienced a higher level of activity in the labor market throughout the pandemic, evidenced by larger values of the baseline transition ratesλ e ,λ u between employed and unemployed groups. Notice also that the estimates for ρ e were positive, whereas the estimates for ρ u were negative, confirming the intuition that an increase in infections has the overall effect of increasing unemployment, in accordance with the model and results of Góes and Gallo (2021). Finally, observe that this sensitivity of unemployment with respect to new infections only begins to be significant towards the start of 2021, as the level of government emergence support in place during most of 2020 resulted in a very weak direct link between infections and unemployment. Figure 4 shows the comparison between the 7-day moving average of nationwide daily numbers of infections, deaths, hospitalizations, ICU admissions, fully vaccinated, as well as estimated unemployed individuals and their corresponding model predictions obtained with calibrated parameters.
To compare the model predictions for the number of unemployed individuals with the reported data, we reverse the steps explained in the previous section to transform the data. Namely, we first obtain a prediction for the monthly proportion of unemployed individuals of working age in each location by using the corresponding daily prediction for the first day of each month. Next, we divide these predictions by the proportion of individuals of working age in each location in 2020 to obtain monthly predictions for the unemployment rate in each location and multiply the result by the size of the labor force for that month, using the size of the labor force in February 2022 for simulations past that date. Finally, we aggregate the resulting number of unemployed individuals for the entire country.
As one can see from Figure 4, all the in-sample model predictions were remarkably adherent to the reports, illustrating the model accuracy when all the available data are used to estimate the time-varying parameters. Specifically, the in-sample calibration error, to be defined in detail in the next section, for the period 19 February 2020 to 28 February 2022 was approximately 11%. In practice, of course, one is interested in out-of-sample predictions, so in the next section we investigate the accuracy of the model when only a portion of the available data are used to estimate the parameters, which are then used to make predictions for subsequent periods.

Backtesting
In order to analyze the forecast capability of the proposed model, we performed a series of backtests. More precisely, we calibrated the model using non-overlapping periods of 20 days each and made out-of-sample predictions 7 days ahead, using the average of the estimated parameter values of the last five days in each period. The comparison between the out-of-sample model predictions and nationwide reports of COVID-19 cases, hospitalizations, ICU admissions, and deaths for some periods can be found in Table 1. Table 1 also presents the calibration and prediction error, which were evaluated as follows: in both cases, we multiplied by 100 the 2 -norm of the difference between the vectors containing the model predictions and reports of daily COVID-19 infections divided by the 2 -norm of the vector containing the reports. The calibration error refers to the error in the in-sample predictions and the prediction error refers to the error in the out-of-sample predictions. Thus, using the notation from Section 3, both the in-sample (Calib. Error) and the out-of-sample (Predic. Error) errors are evaluated using the following formula: We considered only the infections in the evaluation of the errors since we used only the reported COVID-19 cases in the model calibration.   Calibration errors were mainly below 5%, whereas prediction errors were mainly below 20%, which is in agreement with other proposed models, such as . This illustrates the efficiency and the accuracy of the proposed model to design short-term predictions regardless if the incidence is increasing, decreasing, or stable.
We also tested the model prediction capability for longer periods. Figure 5 compares 30-day ahead model predictions with reports during different disease incidence trends. As Figure 5 shows, the model is able to capture the trend in the incidence, extrapolating it for longer periods and providing accurate predictions. The model predictions in the testing set, i.e., the reported data not used in the calibration, are obtained with the average of the estimated parameter values for the last 5 days in the training set, i.e., the dataset used in the calibration. Thus, as long as the observed trend at the end of the training set is kept, the model is able to provide accurate predictions. However, during longer periods, incidence trends can revert, and the model prediction can degenerate as in the predictions for the periods from 19 April to 19 May 2020 and from 1 October to 31 October 2020. Whenever such loss of prediction accuracy occurs, the parameters in the model must be updated. Note that such loss of accuracy in long-term predictions is not unusual due to the nonlinearity and stochastic nature of the phenomena under consideration, as well as the enforcement of new restrictions. For example, since 15 April 2020, international travelers were required to stay in quarantine when entering Canada and since 17 April 2020, air and rail passengers were required to use face masks CIHI (2022); Vogel and Eggertson (2020). Such new rules can be related to the loss of accuracy in the predictions for the period from 19 April to 19 May 2020. In 1 October 2020, stringent containment measures were undertaken in Quebec to control the rise of cases CP24 (2020). Such measures can also be related to the loss of accuracy from 1 October to 31 October 2020.
It is also worth mentioning that we are modeling the COVID-19 outbreak in Canada considering a number of different scales, i.e., different levels of disease severity for each Canadian province or territory, accounting for age range, vaccination, and unemployment. Thus, due to the large complexity of the proposed model and its calibration from real data, we can assert that the performances of the predictions are more than satisfactory.

Long-Term Predictions
It is by now well known that after recovering from the infection or after vaccination, immunity decreases with time. Such loss of immunity can be caused by a series of reasons, including the emergence of new variants of the SARS-CoV-2, such as the delta and omicron variants, which caused large outbreaks worldwide even in countries with efficient vaccination coverage, such as Canada, the United States, and the United Kingdom. Thus, we assumed that after a period of 360 days, anyone who recovered or who was vaccinated returned to the susceptible compartment.
We considered periodic vaccination campaigns with the following periods, 270, 360, 390, 450, 540, and 720 days. Figures 6-8 show the different projections of the evolution of infections, deaths, and unemployed individuals under different vaccination scenarios, using a 10-day average for the time-dependent parameters as of 28 February 2022, namely the end of the calibration period used in Section 4.1.
As can be seen in the graphs, when the vaccination period is shorter than the assumed average time for loss of immunity (namely 360 days), the number of infections rapidly decrease and remain negligible, whereas longer vaccination periods lead to periodic resurgence in infections, with peak numbers comparable to what was observed during the 2020-2022 waves of the pandemic (with the exception of the omicron wave at the start of 2022). Observe that, as the intervals between successive vaccination campaigns increase, the infection waves become less frequent and with larger amplitudes, and similar patterns can be observed in the number of daily deaths. We also observe an increase in the corresponding wave in the monthly unemployment rate as vaccination becomes less frequent, but because of the nonlinear relationship between new infections and the labor market expressed in Equations (26) and (27), their amplitudes first increase and then decrease, with the largest oscillations occurring when the interval between vaccinations is 450 days.

Economic Effect
In this section, we analyze the effect of this pandemic on the economy in a manner similar to Góes and Gallo (2021). Namely, we use the following revised version of Okun's law 10 : where Y is the observed GDP, k is the growth rate of GDP output under constant unemployment, θ is a factor relating changes in unemployment to the growth rate in GDP, and u is the observed unemployment rate. In order to predict the impact of COVID-19 on the economy, we estimated the parameters k, θ using Canadian national data time series. For the unemployment rate u, expressed as a percentage of the labor force, we used quarterly data provided by the OECD 11 and computed the differences between consecutive periods to obtain ∆u. For GDP growth rate ∆Y Y , we used quarterly data provided by the OECD 12 , expressed as percentage change from the previous period.
Applying the ordinary least square (OLS) regression methods, we estimate theses parameters from the first quarter of 2015 to the first quarter of 2022, using data from the preceding 16 quarters for each estimate. The estimated values of the parameters k and θ during the pandemic, together with the corresponding standard deviations, can be found in Table 2, whereas the estimates including the previous 5 years can be found in Figure 9.
As we can see in Table 2, the estimated value for the long-term (i.e., prevailing under constant unemployment) quarterly growth rate in output k dropped from 0.51% to 0.26% from the first to second quarters of 2020 (or from 2.07% to 1.06% in terms of annualized growth rates), illustrating the contraction caused by the onset of the pandemic, and subsequently bounced back to around 0.6% (or a 2.4% annualized growth rate) once the economic stimulus measures were in effect, before stabilizing around 0.4% (or a 1.61% annualized growth rate), which is close to the pre-pandemic 5-year average. Similarly, the estimated value for θ initially increased significantly from 0.71 to 2.74 and then decreased again and stabilized around 1.95, still higher than the pre-pandemic 5-year average of 1.23 but in line with other estimates for Canada and the United States. In other words, the relationship between changes in unemployment and GDP growth first became stronger at the start of the pandemic, before stabilizing towards its previous pattern.  We can now use the results of both the epidemiological model in Equations (7)-(21) and Okun's law to check for consistency between the model and our different data sources, as well as to make long-term predictions of economic effects. Starting with the Canadian unemployment rate, Figure 10 shows the comparison between the quarterly unemployment rate obtained directly from the OECD dataset and the unemployment rate computed from data from Statistics Canada. To compute the Stats Canada unemployment rate, we first aggregate the provincial monthly data used in the previous section into national quarterly averages and then divide the total number of unemployed individuals by the labor force. As we can see, the Stats Canada unemployment rate is considerably noisier than the OECD unemployment rate, primarily because of seasonal adjustment used by the OECD, which is why we chose this time series to estimate Okun's law. Figure 11 shows the comparison between the Stats Canada monthly unemployment rate during the pandemic and the corresponding model predictions for the same period. As expected from the good agreement in the number of unemployed individuals already seen in Figure 4, the unemployment rate predicted by the model is close to the Stats Canada unemployment rate. Accordingly, we can use the model to make the forecasts shown in Figure 11 for the Canadian unemployment rate under different future vaccination scenarios, with a similar pattern as observed for the total number of unemployed individuals in Figure 8.  Moving to economic output, in the top panel of Figure 12 we compare the percentage change in quarterly GDP in Canada obtained directly from the OECD dataset with those obtained using data from Stats Canada by calculating the percentage change in quarterly GDP from monthly GDP data. 13 As we can observe, these two sources provide nearly identical data, and either can be used in comparisons with model predictions. Next, in the bottom panel of Figure 12, we compare the OECD-reported percentage change in quarterly GDP in Canada and those calculated from the estimated Okun law (34) using quarterly changes in the unemployment rate reported by the OECD as an input. We can see that our estimated Okun's law fits the historical data relatively well, including during the pandemic. More specifically, we can observe a precipitous drop in the quarterly percentage change of Canadian GDP in 2020 Q2, coinciding with the beginning of the COVID pandemic in Canada. Subsequently, with the financial support from the Canadian government and control of the first wave in 2020 Q3, the percentage change in GDP experienced a large rebound. Following the second wave of transmission, the speed of economic recovery began to slow down again in 2020 Q4 and stabilized at historical levels from 2021 Q1 onward. Figure 13 shows the predictions of percentage change in quarterly GDP in Canada obtained from the estimated 14 Okun's law (34) using the quarterly changes in unemployment rate predicted by the model as an input. As we can see, the quarterly changes in GDP predicted in this way also fit the data well during the pandemic. Accordingly, we use the model to make forecasts for quarterly growth rates in GDP under different future vaccination scenarios. Unsurprisingly, given the pattern of unemployment observed in Figure 11, we also see a similar behavior for the cycles of GDP growth, namely less frequent vaccination campaigns lead to GDP growth oscillations of increasing periods, but increasing and then decreasing amplitudes.

Conclusions and Further Research Directions
In this paper, we developed a combined economic-epidemiological model to investigate the effects of the COVID-19 pandemic in Canada. The epidemiological part of the model consists of an extension of the susceptible-exposed-infected-removed (SEIR) along the lines of ,i.e., incorporating time-varying parameters (when applicable), multiple geographic regions (Canadian provinces and territories), and multiple age groups. This last feature was introduced in the model so that we could identify a working age group and establish a link between the pandemic and the economy via unemployment. For this, we further subdivided each compartment in the model into the employed and unemployed sub-compartments and adjusted the epidemiological dynamics accordingly, for example, by assuming that contact rates between individuals in unemployed compartments are lower than in employed ones, reflecting less mobility between home and work as well as less exposure to the general public for those who are unemployed. Crucially, we also assumed that the labor market dynamics were also affected by the pandemic through time-varying transition rates between employment status with an explicit dependence on new infections, along the lines suggested in Góes and Gallo (2021). To explore the economic effects further, we assumed an Okun's law relationship between the changes in unemployment and GDP, and estimated them before and during the pandemic.
As shown in Figure 4, the calibrated model was able to fit the observed data remarkably well, with an in-sample calibration error of approximately 11% when the entire available dataset was used in the estimation of the parameters. The performance was also satisfactory when shorter periods (20 days) were used in the estimation and we used the model to make predictions 7 days ahead, with calibration errors mainly below 5% and prediction errors mainly below 20%, as shown in Table 1. However, Figure 5 illustrates that, as should be expected, the performance of the model deteriorates when it is used to make longer term (i.e., 30 days ahead) predictions during periods of rapid changes in behavior, and consequently model parameters, such as when new lockdown measures were imposed. Even in these instances, however, periodic re-calibration of the underlying parameters is enough to restore the performance of the model.
With this caveat in mind, we used the model to make long-term predictions (i.e., up to 2030) for COVID-19 under different future scenarios for vaccinations and assuming that loss of immunity occurs on average 360 days after either vaccine or recovery from the disease. Figures 6-8 show the predictions for infections, deaths, and unemployed individuals for a variety of scenarios, as an example of how the model can be used by policymakers interested in the design of vaccination campaigns.
Regarding the economic impact, Figure 9 shows that estimates for long-term quarterly growth rate at constant unemployment, namely the parameter k in Okun's law (34), exhib-ited a sharp decline at the start of the pandemic, followed by a quick rebound throughout 2020 and a reversion to the pre-pandemic average toward the end of 2021. On the other hand, the response of quarterly growth rates to unemployment, namely the parameter θ in Okun's law (34), became stronger at the start of the pandemic and then stabilized at a level higher than the pre-pandemic average. Interestingly, whereas the uncertainty in the estimate of k increased since the start of the pandemic, the uncertainty in the estimate of θ appears to have decreased.
The links between the labor market and disease dynamics embedded throughout the model, most notably in Equations (26) and (27), allowed us to make the long-term forecasts for unemployment rates, shown in Figure 11, with cycles following a pattern differently from that observed for infections and deaths, namely as the intervals between vaccination campaigns increased, the amplitude of oscillations in unemployment rates first became larger, but eventually began to decrease. This in turn can be converted into long-term forecasts for the GDP growth rates shown in Figure 13.
As mentioned above, we offer these long-term predictions as illustrations of how the model can be used as a tool to aid policy-making, rather than authoritative statements about the future dynamics of these variables. Specifically, what we have in mind is that policy makers should run the model under a variety of plausible parameter values and policy scenarios (e.g., vaccination intervals) to obtain the range of outcomes for the variables of interest, notably deaths and unemployed individuals.
Having explored here a reasonably detailed epidemiological model for Canada, including multiple compartments for different stages of the disease (e.g., moderate, hospitalized, admitted to ICU), as well as age groups and geographic locations, in future work we plan to focus on simpler versions of the model but extend the analysis to a larger number of countries. Specifically, we will consider aggregate national data and a simplified epidemiological structure for a set of representative countries (e.g., USA, United Kingdom, Italy, China, etc.) and explore the differences in the infection-unemployment nexus in the context of a reduced model, which opens up the possibility of analytical results, in addition to the estimation and simulation-based results presented here.
Author Contributions: All the authors contributed equally to the development, the model, and the interpretation of the results. V.A. conducted the estimation and simulation of the epidemiological part of the model and W.P. conducted the estimation and simulation of the economic part of the model. All authors have read and agreed to the published version of the manuscript.
Funding: V.A. acknowledges the financial support from Fundação Butantan and Fundação de Amparo à Pesquisa e Inovação do Estado de Santa Catarina through grant nos. 01/2020 and 00002847/2021, respectively. J.Z. acknowledges the financial support from Khalifa University, CNPq, and Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro through grant nos. FSU-2020-09, 307873/2013-7, andE-26/202.927/2017, respectively. M.G. and W.P. acknowledge the financial support from the Natural Sciences and Engineering Research Council of Canada through the Discovery Grants program.

Data Availability Statement:
Publicly available datasets were used in this study. All code and data needed to implement the model can be found at https://github.com/viniciusalbani/Covideconomics (accessed on 5 October 2022).

Conflicts of Interest:
The authors declare no competing interests.
2 Data on air travel distance downloaded 20 November 2020 from https://www.distancefromto.net/ (accessed on 20 November 2020). 3 For clarity, consider provinces 1, 2, 3 and suppose that the distance between the centers of 1 and 2 is in the same range as the distance between the centers of 2 and 3, say 3000 to 4000 K. We then impose the constraint b 1,2 = b 2,3 = b 4 and estimate b 4 using data for both pairs of provinces. This reduces the original 78 parameters b i,j , i = j, to the five parameters b 1 , . . . , b 5 . 4 The term 1 inside the argument for the logarithms is added to avoid indefinite values when the number of reported infections is zero. 5 Notice that, to find the values for j = 1,2,3, we just have to consider the series values for previous dates. 6 Data downloaded on 15 March 2022 from https://resources-covid19canada.hub.arcgis.com/datasets/provincial-daily-totals/ explore (accessed on 15 March 2022). 7 Data downloaded on 20 Novemver 2020 from https://www.statcan.gc.ca/eng/subjects-start/population_and_demography (accessed on 20 Novemver 2020). 8 Data downloaded on 1 July 2022 from https://www150.statcan.gc.ca/t1/tbl1/en/tv.action?pid=1410028703 (accessed on 1 July 2022). 9 In accordance with data in https://www.cdc.gov/coronavirus/2019-ncov/variants/omicron-variant.html (accessed on 20 May 2022). 10 The original version proposed in Okun (1963) has the causal relationship running from the growth rate of a national gross domestic product (GDP) to changes in unemployment rate.