Covid-19 SIR Modeling Malaysia

Description

This SIR model is modeled by finding the best parameters that fit to Covid-19 Malaysian data up to Fri, 30 Jul 2021.

The SIR model is defined by the following set of ordinary differential equations (Ref: Wikipedia):

\[ \frac{dS}{dt} = -\frac{\beta IS}{N} \\ \frac{dI}{dt} = \frac{\beta IS}{N} - \gamma I \\ \frac{dR}{dt} = \gamma I \] where \(S\) = Susceptible population, \(I\) = Infected (active), \(R\) = Recovered (recover + death), and \(N = S(t) + I(t) + R(t)\). \(\gamma\) is defined as \(\frac{1}{T_r}\), where \(T_r\) = time until recovery, which can be obtained from clinical observation. \(\beta\) is defines as \(\frac{1}{T_c}\), where \(T_c\) = time between contacts.

The basic reproduction number, \(R_0\) is defined as:

\[ R_0 = \frac{\beta}{\gamma} \] \(R_0\) can also be seen as \(\frac{T_r}{T_c}\), i.e. the number of contacts before someone recovers, during which he is infectious.

Apart from finding the \(\beta\) and \(\gamma\) parameter values, experts are concerned about the number of susceptible individuals. Thus, the model is also iterated over decreasing number of initial susceptible individuals (\(n\)) to find this optimal \(n\).

To fit the model to the data, the following weighted loss function is defined, extending the basic residual sum of squares (RSS) formulation:

\[ RSS = w_1\sum_{i=1}^t(I_i - \hat I_i)^2 + w_2\sum_{i=1}^t(R_i - \hat R_i)^2 + w_3\sum_{i=1}^t([I_i + R_i] - [\hat I_i + \hat R_i])^2 \] where \(w_1\), \(w_2\) and \(w_3\) allow adjustment of weightage given to \(I\), \(R\) respectively to balance the model fit. We are inclined to give more weight to \(I\) through \(w_1\) and \(I + R\) through \(w_3\).

Here, given the limitation of the optimization, one might argue that the \(S\) found might instead be \(E\) in SEIR model. This could be explored further in SEIR model, by allowing incubation period in the model. This is beyond the present SIR model to optimize.

Libraries

library(tidyverse)
library(deSolve)
library(lubridate)
library(ggpubr)
library(plotly)

Data

Sourced from local datasets.

# mys_data = read_csv("https://wnarifin.github.io/covid-19-malaysia/covid-19_my_full.csv") %>% select(1:11) %>% rename(Date = "date")
mys_data = read_csv("covid-19_my_full.csv") %>% select(1:11)  %>% rename(Date = "date")  # comment out to use local data
# rename date to Date for compatibility with original codes

# add active variable
mys_data$active = with(mys_data, total_cases - total_deaths - total_recover)
names(mys_data)

##  [1] "Date"           "location"       "new_cases"      "new_deaths"    
##  [5] "total_cases"    "total_deaths"   "recover"        "total_recover" 
##  [9] "icu"            "support"        "imported_cases" "active"

Functions

# SIR model
# S - Susceptible, I - Infected, R - Removed (recovered + death), n - initial susceptible population
SIR = function(time, state, parameters) {
  par = as.list(c(state, parameters))
  with(par, {
    dS = (-beta * I * S) / n
    dI = ((beta * I * S) / n) - (gamma * I)
    dR = gamma * I
    list(c(dS, dI, dR))
  })
}

# RSS 
RSS = function(parameters) {
  names(parameters) = c("beta","gamma")
  out = ode(y=init, times=Day, func=SIR, parms=parameters)
  fit1 = out[,3]
  fit2 = out[,4]
  # give weight to less well predicted curve
  w1 = 3
  w2 = 1
  w3 = 3
  # RSS Original scale
  # w1*sum((Active - fit1)^2) + w2*sum((Removed - fit2)^2) # find min RSS on original scale, weighted. Works without too much strict setting for gamma.
  w1*sum((Active - fit1)^2) + w2*sum((Removed - fit2)^2) + w3*sum(((Active+Removed) - (fit1+fit2))^2)
  # ---
  # RSS log scale
  # w1*sum((log(Active) - log(fit1))^2) + w2*sum((log(Removed) - log(fit2))^2)  # find min RSS on log scale, weighted. Needs more work on correct gamma.
  # w1*sum((log(Active) - log(fit1))^2) + w2*sum((log(Removed) - log(fit2))^2) + w3*sum((log(Active+Removed) - log(fit1+fit2))^2)
  # Give balance to prediction of small values and large values.
  # ---
  # MAPE
  # w1*(mean(abs(log(Active) - log(fit1))/log(Active)))  + w2*(mean(abs(log(Removed) - log(fit2))/log(Removed)))  # MAPE
  # w1*(mean(abs(Active - fit1)/Active) / length(fit1)) + w2*(sum(abs(Removed - fit2)/Removed)) + w3*(mean(abs((Active+Removed) - (fit1+fit2))/(Active+Removed)))
  # Intuitive, but may be problematic when denominator = 0, or error > 1 for large prediction.
}

Parameters

# Uncomment under each header to perform SIR for the period

# Pre MCO #
# name = "Pre MCO"
# start_date  = "2020-03-01"
# end_date    = "2020-03-17"

# MCO All -> Today #
# name = "MCO All"
# start_date  = "2020-03-18"
# end_date    = today()  # will analyze up to max available date

# MCO week 2 -> End MCO 4 # to take into account MCO effect after 1 week
# name = "MCO week 2 -> End MCO 4"
# start_date  = "2020-03-25"
# end_date    = "2020-05-03"

# MCO 1 #
# name = "MCO 1"
# start_date  = "2020-03-18"
# end_date    = "2020-03-31"

# MCO 2 #
# name = "MCO 2"
# start_date  = "2020-04-01"
# end_date    = "2020-04-14"

# MCO 3 #
# name = "MCO 3"
# start_date  = "2020-04-15"
# end_date    = "2020-04-28"

# MCO 4 #
# name = "MCO 4"
# start_date  = "2020-04-28"
# end_date    = "2020-05-03"

# CMCO #
# name = "CMCO"
# start_date  = "2020-05-04"
# end_date    = "2020-06-09"

# RMCO 1 #
# name = "RMCO 1"
# start_date  = "2020-06-10"
# end_date    = "2020-08-31"

# RMCO 2 #
# name = "RMCO 2"
# start_date  = "2020-09-01"
# end_date    = today()  # will analyze up to max available date

# Sabah Election #
# name = "Sabah Election"
# start_date  = "2020-09-26"
# end_date    = today()  # will analyze up to max available date

# Sabah Election -> MCO 2.0 # See R0 w/out MCO
# name = "Sabah Election"
# start_date  = "2020-09-26"
# end_date    = "2021-01-12"  # up to MCO 2

# Sabah Election -> Unrestricted Travel # See R0 until 2020-12-06
# name = "Sabah Election"
# start_date  = "2020-09-26"
# end_date    = "2020-12-06"  # up to MCO 2

# RMCO 2 Unrestricted Travel -> Today #
# name = "RMCO 2 Unrestricted Travel"
# start_date  = "2020-12-07"
# end_date    = today()  # will analyze up to max available date

# RMCO 2 Unrestricted Travel -> MCO 2.0 # See R0 w/out MCO
# name = "RMCO 2 Unrestricted Travel"
# start_date  = "2020-12-07"
# end_date    = "2021-01-12"  # up to MCO 2

# MCO 2.0 #
# name = "MCO 2"
# start_date  = "2021-01-13"
# end_date    = today()  # will analyze up to max available date

# Opening of school
# Pre-school-std 2 1-3-21
# Std 3 - 6 8-3-21
# Secondary school 4-4-21
name = "Secondary School Reopened"
# start_date  = "2021-03-01"  # even too early, active cases going down
# start_date  = "2021-03-08"  # too early, active cases going down
start_date  = "2021-04-04"  # notable increase in cases, esp in Kelantan
end_date    = today()  # will analyze up to max available date

# Basic info from data
Infected   = mys_data %>% filter(Date>=ymd(start_date), Date<=ymd(end_date)) %>% pull(total_cases)
Recovered  = mys_data %>% filter(Date>=ymd(start_date), Date<=ymd(end_date)) %>% pull(total_recover)
Death      = mys_data %>% filter(Date>=ymd(start_date), Date<=ymd(end_date)) %>% pull(total_deaths)
Active     = Infected - Recovered - Death
Removed    = Recovered + Death
Date       = mys_data %>% filter(Date>=ymd(start_date), Date<=ymd(end_date)) %>% pull(Date)
Day        = 1:(length(Date))

Run Optimization for SIR model

# SIR parameters
# parameters_values = c(1/2, 1/14)  # set reasonable start, R0 = 7
# parameters_values_lower = c(1/100, 1/42)  # days recover max 6 weeks
# parameters_values_lower = c(1/100, 1/18.5)  # setting for log scale, days recover 6 weeks, reduce bound bcs not many severe, observed show quick recovery
# parameters_values_upper = c(1, 1/11)  # updated to min 11 days
# based on previous models in paper 1 betas (0.273, 0.152) gammas (0.033, 0.053)
parameters_values = c(0.273, 0.033)
parameters_values_lower = c(0.1, 0.01)
# parameters_values_upper = c(0.35, 0.07)  # max gamma 0.07 i.e. 14 days
parameters_values_upper = c(1, 0.07)  # max gamma 0.07 i.e. 14 days
# parameters_values_upper = c(0.5, 0.1)  # allows freedom in parameter space exploration, gamma 1/10
# cannot allow this bcs it will never converge as N cycles between 32mil (odd max step) & 64mil (even max step)!

# Placeholder to find optimal susceptible population
max_step = 7   # 7 seems sufficient
N_in_steps = data.frame(step=1:max_step, N=rep(NA,max_step), Loc=rep(NA,max_step))

# Initial values
N           = 32.68E6/4  # 4th quarter, 2019
Opt_min_loc = 1  # Optimum minimum n location in output vector
step_i      = 1  # initialize step

# Steps
for (step_i in 1:max_step) {
  cat("=== Step ", step_i, ": Finding optimal values of beta, gamma and n ===\n", sep ="")
  N = N[Opt_min_loc]  # Susceptible population from previous Step
  if (N > 32.68E6) {N = 32.68E6/4}  # avoid N > Malaysian Population
  # initial susceptible population of Malaysia
  p_start = 0.05
  p_end = 4  # max p of population, also include up direction
  p_step = 0.05
  susceptible = seq(p_start, p_end, p_step)
  N = N * susceptible
  inits = data.frame(S=N-Infected[1]-Recovered[1]-Death[1], I=Infected[1]-Death[1]-Recovered[1], R=Recovered[1]+Death[1])
  Opts = vector("list", length(N))
  for (i in 1:length(N)) {
    n = N[i]
    init = setNames(as.numeric(inits[i,]), c("S", "I", "R"))
    Opt_ = optim(parameters_values, RSS, method = "L-BFGS-B", lower = parameters_values_lower, upper = parameters_values_upper)
    Opts[i]  = list(Opt_)
  }
  Opt_value = sapply(Opts, function(x) x$value)
  Opt_min_loc = which(Opt_value == min(Opt_value))
  N_in_steps[step_i, "N"] = N[Opt_min_loc]
  N_in_steps[step_i, "Loc"] = Opt_min_loc
  cat("=== Found n = ", N[Opt_min_loc], " at location ", Opt_min_loc, " in vector N ===\n\n", sep = "")
  step_i = step_i + 1
  if (step_i == max_step + 1) {
    cat("=== Finalizing results =====\n")
    cat("============================\n")
    print(N_in_steps)
  }
}

## === Step 1: Finding optimal values of beta, gamma and n ===
## === Found n = 2859500 at location 7 in vector N ===
## 
## === Step 2: Finding optimal values of beta, gamma and n ===
## === Found n = 2716525 at location 19 in vector N ===
## 
## === Step 3: Finding optimal values of beta, gamma and n ===
## === Found n = 2716525 at location 20 in vector N ===
## 
## === Step 4: Finding optimal values of beta, gamma and n ===
## === Found n = 2716525 at location 20 in vector N ===
## 
## === Step 5: Finding optimal values of beta, gamma and n ===
## === Found n = 2716525 at location 20 in vector N ===
## 
## === Step 6: Finding optimal values of beta, gamma and n ===
## === Found n = 2716525 at location 20 in vector N ===
## 
## === Step 7: Finding optimal values of beta, gamma and n ===
## === Found n = 2716525 at location 20 in vector N ===
## 
## === Finalizing results =====
## ============================
##   step       N Loc
## 1    1 2859500   7
## 2    2 2716525  19
## 3    3 2716525  20
## 4    4 2716525  20
## 5    5 2716525  20
## 6    6 2716525  20
## 7    7 2716525  20

# Saving optimized parameters
Opt = Opts[[Opt_min_loc]]
Opt$message  # make sure converge

## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Opt_par = setNames(Opt$par, c("beta", "gamma"))
R0 = (Opt_par['beta']/Opt_par['gamma']); names(R0) = "R0"
parameters_values_lower; parameters_values; parameters_values_upper  # just to check whether values in range

## [1] 0.10 0.01

## [1] 0.273 0.033

## [1] 1.00 0.07

# Print final parameters
cat("beta = ", Opt_par[['beta']], ", infectious contact rate (/person/day)\n",
    "gamma = ", Opt_par[['gamma']], ", recovery rate (/day)\n",
    "R_0 = ", R0, " number infected/person\n",
    "Recovery days = ", 1/Opt_par[['gamma']], " days",
    sep = "")

## beta = 0.1346302, infectious contact rate (/person/day)
## gamma = 0.07, recovery rate (/day)
## R_0 = 1.923289 number infected/person
## Recovery days = 14.28571 days

Fit Data

# time in days for fitting
t = 1:max(Day)
n = N[Opt_min_loc]
init = setNames(as.numeric(inits[Opt_min_loc,]), c("S", "I", "R"))
# get the fitted values from our SIR model
fitted_projected = data.frame(ode(y=init, times=t, func=SIR, parms=Opt_par))
# add Date, Active, Removed
fitted_projected$Date = Date
fitted_projected$A = Active
fitted_projected$Rm = Removed
fitted_projected

##     time       S         I        R       Date      A     Rm
## 1      1 2029116  14509.00 336450.0 2021-04-04  14509 336450
## 2      2 2027635  14958.74 337481.3 2021-04-05  14278 337751
## 3      3 2026109  15421.27 338544.6 2021-04-06  14161 339168
## 4      4 2024538  15896.86 339640.6 2021-04-07  14097 340371
## 5      5 2022919  16385.84 340770.4 2021-04-08  14203 341550
## 6      6 2021252  16888.48 341934.9 2021-04-09  14805 342802
## 7      7 2019535  17405.08 343135.1 2021-04-10  15059 344058
## 8      8 2017767  17935.93 344372.0 2021-04-11  15574 345282
## 9      9 2015947  18481.33 345646.5 2021-04-12  15835 346338
## 10    10 2014074  19041.57 346959.7 2021-04-13  16300 347640
## 11    11 2012145  19616.95 348312.7 2021-04-14  16696 349133
## 12    12 2010161  20207.75 349706.5 2021-04-15  17575 350402
## 13    13 2008119  20814.28 351142.1 2021-04-16  18600 351928
## 14    14 2006017  21436.80 352620.8 2021-04-17  19094 353765
## 15    15 2003856  22075.61 354143.7 2021-04-18  19854 355200
## 16    16 2001632  22730.99 355711.8 2021-04-19  20522 356610
## 17    17 1999345  23403.21 357326.4 2021-04-20  21268 358205
## 18    18 1996994  24092.54 358988.6 2021-04-21  21687 360126
## 19    19 1994576  24799.24 360699.8 2021-04-22  22014 362674
## 20    20 1992090  25523.57 362461.0 2021-04-23  22512 365023
## 21    21 1989536  26265.77 364273.5 2021-04-24  22926 367326
## 22    22 1986910  27026.09 366138.6 2021-04-25  23753 369189
## 23    23 1984213  27804.75 368057.6 2021-04-26  24713 371005
## 24    24 1981441  28601.97 370031.7 2021-04-27  25414 373037
## 25    25 1978595  29417.95 372062.3 2021-04-28  26719 374874
## 26    26 1975671  30252.89 374150.6 2021-04-29  28093 376832
## 27    27 1972670  31106.95 376298.1 2021-04-30  29227 379486
## 28    28 1969589  31980.31 378506.1 2021-05-01  29631 381963
## 29    29 1966426  32873.09 380775.8 2021-05-02  30339 384673
## 30    30 1963181  33785.44 383108.7 2021-05-03  30753 386759
## 31    31 1959851  34717.45 385506.2 2021-05-04  31516 389116
## 32    32 1956436  35669.21 387969.7 2021-05-05  32939 391437
## 33    33 1952934  36640.79 390500.4 2021-05-06  33762 394165
## 34    34 1949343  37632.21 393099.8 2021-05-07  34789 397636
## 35    35 1945662  38643.49 395769.4 2021-05-08  36564 400380
## 36    36 1941890  39674.62 398510.4 2021-05-09  37060 403617
## 37    37 1938025  40725.56 401324.3 2021-05-10  37396 407088
## 38    38 1934066  41796.23 404212.4 2021-05-11  38499 409958
## 39    39 1930012  42886.54 407176.2 2021-05-12  40101 413121
## 40    40 1925862  43996.34 410217.0 2021-05-13  41582 416495
## 41    41 1921613  45125.46 413336.1 2021-05-14  41471 420719
## 42    42 1917266  46273.71 416535.0 2021-05-15  42135 424195
## 43    43 1912819  47440.84 419814.9 2021-05-16  41889 428221
## 44    44 1908271  48626.56 423177.1 2021-05-17  43506 431050
## 45    45 1903621  49830.56 426623.0 2021-05-18  44827 434594
## 46    46 1898869  51052.48 430153.8 2021-05-19  47340 438156
## 47    47 1894012  52291.91 433770.8 2021-05-20  50171 442131
## 48    48 1889051  53548.40 437475.1 2021-05-21  52106 446689
## 49    49 1883986  54821.46 441268.0 2021-05-22  53682 451433
## 50    50 1878814  56110.54 445150.5 2021-05-23  57022 455069
## 51    51 1873536  57415.08 449123.8 2021-05-24  60018 458582
## 52    52 1868152  58734.42 453188.9 2021-05-25  63458 462431
## 53    53 1862660  60067.88 457346.9 2021-05-26  66208 467159
## 54    54 1857061  61414.75 461598.8 2021-05-27  69408 471816
## 55    55 1851355  62774.22 465945.3 2021-05-28  72823 476691
## 56    56 1845542  64145.47 470387.4 2021-05-29  76218 482316
## 57    57 1839621  65527.63 474925.9 2021-05-30  78017 487516
## 58    58 1833594  66919.74 479561.5 2021-05-31  79523 492834
## 59    59 1827459  68320.84 484294.9 2021-06-01  80474 498988
## 60    60 1821218  69729.89 489126.6 2021-06-02  82274 504891
## 61    61 1814872  71145.80 494057.2 2021-06-03  83331 512043
## 62    62 1808420  72567.45 499087.2 2021-06-04  84369 518753
## 63    63 1801865  73993.65 504216.8 2021-06-05  85607 524967
## 64    64 1795205  75423.17 509446.4 2021-06-06  86628 530187
## 65    65 1788444  76854.75 514776.1 2021-06-07  84269 537817
## 66    66 1781582  78287.07 520206.0 2021-06-08  82797 544855
## 67    67 1774620  79718.77 525736.3 2021-06-09  81575 552316
## 68    68 1767560  81148.45 531366.6 2021-06-10  79848 559714
## 69    69 1760403  82574.68 537097.0 2021-06-11  78864 567547
## 70    70 1753152  83995.98 542927.0 2021-06-12  76247 575957
## 71    71 1745808  85410.84 548856.3 2021-06-13  73324 584184
## 72    72 1738373  86817.74 554884.3 2021-06-14  71625 590832
## 73    73 1730849  88215.10 561010.5 2021-06-15  70112 597764
## 74    74 1723239  89601.35 567234.2 2021-06-16  67949 605077
## 75    75 1715546  90974.87 573554.4 2021-06-17  66097 612667
## 76    76 1707771  92334.05 579970.3 2021-06-18  65602 619602
## 77    77 1699917  93677.25 586480.8 2021-06-19  64523 626592
## 78    78 1691987  95002.82 593084.7 2021-06-20  63815 632593
## 79    79 1683985  96309.12 599780.7 2021-06-21  62918 638101
## 80    80 1675913  97594.50 606567.5 2021-06-22  62027 643735
## 81    81 1667774  98857.31 613443.4 2021-06-23  60816 650190
## 82    82 1659572 100095.92 620407.0 2021-06-24  61162 655685
## 83    83 1651310 101308.71 627456.3 2021-06-25  60117 662542
## 84    84 1642991 102494.07 634589.5 2021-06-26  60646 667816
## 85    85 1634620 103650.42 641804.8 2021-06-27  61395 672653
## 86    86 1626199 104776.21 649099.9 2021-06-28  61812 677454
## 87    87 1617732 105869.91 656472.7 2021-06-29  62844 682859
## 88    88 1609224 106930.04 663920.9 2021-06-30  64129 687850
## 89    89 1600678 107955.15 671442.1 2021-07-01  65453 693514
## 90    90 1592097 108943.85 679033.8 2021-07-02  66084 699865
## 91    91 1583487 109894.77 686693.3 2021-07-03  66958 705649
## 92    92 1574850 110806.62 694418.1 2021-07-04  67669 710983
## 93    93 1566192 111678.16 702205.3 2021-07-05  69447 715592
## 94    94 1557515 112508.22 710052.1 2021-07-06  72201 720492
## 95    95 1548824 113295.67 717955.5 2021-07-07  74344 725446
## 96    96 1540123 114039.48 725912.5 2021-07-08  77275 731383
## 97    97 1531416 114738.66 733920.0 2021-07-09  80665 737173
## 98    98 1522708 115392.32 741974.8 2021-07-10  84021 743170
## 99    99 1514002 115999.63 750073.8 2021-07-11  87841 748455
## 100  100 1505301 116559.85 758213.7 2021-07-12  91272 753598
## 101  101 1496612 117072.32 766391.1 2021-07-13  96236 759713
## 102  102 1487936 117536.45 774602.7 2021-07-14 101359 766208
## 103  103 1479278 117951.74 782845.0 2021-07-15 108369 772413
## 104  104 1470642 118317.78 791114.8 2021-07-16 114053 779270
## 105  105 1462032 118634.25 799408.4 2021-07-17 119814 786037
## 106  106 1453452 118900.91 807722.4 2021-07-18 124593 791968
## 107  107 1444904 119117.59 816053.3 2021-07-19 128997 798536
## 108  108 1436393 119284.23 824397.7 2021-07-20 133703 806196
## 109  109 1427922 119400.85 832751.9 2021-07-21 137587 814297
## 110  110 1419495 119467.54 841112.6 2021-07-22 142051 822867
## 111  111 1411114 119484.50 849476.2 2021-07-23 147386 833105
## 112  112 1402784 119451.98 857839.3 2021-07-24 153633 842760
## 113  113 1394506 119370.33 866198.4 2021-07-25 160903 852535
## 114  114 1386285 119239.97 874550.0 2021-07-26 165840 862114
## 115  115 1378123 119061.42 882890.8 2021-07-27 170224 873847
## 116  116 1370022 118835.22 891217.5 2021-07-28 175113 886363
## 117  117 1361986 118562.04 899526.7 2021-07-29 179179 899467

Fit measure, original scales

# log scale
# tss1 = sum((log(fitted_projected$A) - mean(log(fitted_projected$A)))^2); tss1
# rss1 = sum((log(fitted_projected$A) - log(fitted_projected$I))^2); rss1
# R2_1 = 1 - (rss1 / tss1); R2_1
# tss2 = sum((log(fitted_projected$Rm) - mean(log(fitted_projected$Rm)))^2); tss2
# rss2 = sum((log(fitted_projected$Rm) - log(fitted_projected$R))^2); rss2
# R2_2 = 1 - (rss2 / tss2); R2_2
# mape1 = sum(abs(log(fitted_projected$A) - log(fitted_projected$I))/log(fitted_projected$A)) / length(fitted_projected$A); mape1
# mape2 = sum(abs(log(fitted_projected$Rm) - log(fitted_projected$R))/log(fitted_projected$Rm)) / length(fitted_projected$Rm); mape2
# original scale
tss1 = sum((fitted_projected$A - mean(fitted_projected$A))^2); tss1

## [1] 181269523113

rss1 = sum((fitted_projected$A - fitted_projected$I)^2); rss1

## [1] 57022729369

R2_1 = 1 - (rss1 / tss1); R2_1

## [1] 0.6854257

tss2 = sum((fitted_projected$Rm - mean(fitted_projected$Rm))^2); tss2

## [1] 3.307845e+12

rss2 = sum((fitted_projected$Rm - fitted_projected$R)^2); rss2

## [1] 36179418452

R2_2 = 1 - (rss2 / tss2); R2_2

## [1] 0.9890625

rmse1 = sqrt(mean((fitted_projected$A - fitted_projected$I)^2)); rmse1

## [1] 22076.54

rmse2 = sqrt(mean((fitted_projected$Rm - fitted_projected$R)^2)); rmse2

## [1] 17584.82

mape1 = mean(abs(fitted_projected$A - fitted_projected$I)/fitted_projected$A); mape1

## [1] 0.2168597

mape2 = mean(abs(fitted_projected$Rm - fitted_projected$R)/fitted_projected$Rm); mape2

## [1] 0.02198104

Plots

# color settings
colors = c("Susceptible" = "black", "Recovered" = "green", "Infectious" = "red", 
           "Observed Active" = "orange", "Observed Recovered" = "blue")

# plot fit the data
fitplot1 = ggplot(fitted_projected, aes(x = Date)) + geom_line(aes(y = I, color = "Infectious")) + 
  geom_point(aes(y = A, color = "Observed Active")) + geom_line(aes(y = R, color = "Recovered")) + 
  geom_point(aes(y = Rm, color = "Observed Recovered")) +
  labs(y = "Number of cases", title = paste("COVID-19 SIR model Malaysia, fitted and observed", name),
       subtitle = paste("Projection from data:", start_date, "to", fitted_projected$Date[max(Day)]),
       color = paste0("Model fit:\n",
                      "R square (I) = ", round(R2_1,3), "\n",
                      "R square (R) = ", round(R2_2,3), "\n",
                      "MAPE (I) = ", round(mape1,3), "\n",
                      "MAPE (R) = ", round(mape2,3), "\n\n",
                      "SIR parameters:\n",
                      "R0 = ", round(R0, 3), "\n",
                      "beta = ", round(Opt_par[1], 3), "\n",
                      "gamma = ", round(Opt_par[2], 3), "\n\n",
                      "Susceptible = ", round(n), "\n",
                      "Peak Active = ", round(max(fitted_projected$I)), "\n")) +
  scale_colour_manual(values = colors)
fitplot1

# uncomment ggsave line to save plot in png format, make sure folder "plots_srn" is created
# ggsave(paste0("plots_srn/fit", name, ".png"), width = 12, height = 9)

# plot fit the data, in log10
fitplot1_log = ggplot(fitted_projected, aes(x = Date)) + geom_line(aes(y = I, color = "Infectious")) + 
  geom_point(aes(y = A, color = "Observed Active")) + geom_line(aes(y = R, color = "Recovered")) + 
  geom_point(aes(y = Rm, color = "Observed Recovered")) + scale_y_log10() + 
  labs(y = "Number of cases", title = paste("COVID-19 SIR model Malaysia, fitted and observed,", name, "log10"),
       subtitle = paste("Projection from data:", start_date, "to", fitted_projected$Date[max(Day)]),
       color = paste0("Model fit:\n",
                      "R square (I) = ", round(R2_1,3), "\n",
                      "R square (R) = ", round(R2_2,3), "\n",
                      "MAPE (I) = ", round(mape1,3), "\n",
                      "MAPE (R) = ", round(mape2,3), "\n\n",
                      "SIR parameters:\n",
                      "R0 = ", round(R0, 3), "\n",
                      "beta = ", round(Opt_par[1], 3), "\n",
                      "gamma = ", round(Opt_par[2], 3), "\n\n",
                      "Susceptible = ", round(n), "\n",
                      "Peak Active = ", round(max(fitted_projected$I)), "\n")) +
  scale_colour_manual(values = colors)  
fitplot1_log

# uncomment ggsave line to save plot in png format, make sure folder "plots_srn" is created
# ggsave(paste0("plots_srn/fit", name, "_log.png"), width = 12, height = 9)

Projected Data

# last date to project
last_date = "2022-06-30"
# time in days for predictions
t = 1:as.integer(ymd(last_date) + 1 - ymd(start_date))
# get the fitted values from our SIR model
fitted_projected = data.frame(ode(y=init, times=t, func=SIR, parms=Opt_par))
# add add Date, Active, Removed
fitted_projected$Date = ymd(start_date) + days(t - 1)
fitted_projected$A = c(Active, rep(NA, length(t) - length(Active)))
fitted_projected$Rm = c(Removed, rep(NA, length(t) - length(Active)))
head(fitted_projected, 10); tail(fitted_projected, 10)

##    time       S        I        R       Date     A     Rm
## 1     1 2029116 14509.00 336450.0 2021-04-04 14509 336450
## 2     2 2027635 14958.74 337481.3 2021-04-05 14278 337751
## 3     3 2026109 15421.27 338544.6 2021-04-06 14161 339168
## 4     4 2024538 15896.86 339640.6 2021-04-07 14097 340371
## 5     5 2022919 16385.84 340770.4 2021-04-08 14203 341550
## 6     6 2021252 16888.48 341934.9 2021-04-09 14805 342802
## 7     7 2019535 17405.08 343135.1 2021-04-10 15059 344058
## 8     8 2017767 17935.93 344372.0 2021-04-11 15574 345282
## 9     9 2015947 18481.33 345646.5 2021-04-12 15835 346338
## 10   10 2014074 19041.57 346959.7 2021-04-13 16300 347640

##     time        S        I       R       Date  A Rm
## 444  444 908408.9 95.61234 1471570 2022-06-21 NA NA
## 445  445 908404.7 93.25356 1471577 2022-06-22 NA NA
## 446  446 908400.5 90.95296 1471584 2022-06-23 NA NA
## 447  447 908396.5 88.70909 1471590 2022-06-24 NA NA
## 448  448 908392.5 86.52057 1471596 2022-06-25 NA NA
## 449  449 908388.7 84.38602 1471602 2022-06-26 NA NA
## 450  450 908385.0 82.30412 1471608 2022-06-27 NA NA
## 451  451 908381.3 80.27357 1471613 2022-06-28 NA NA
## 452  452 908377.7 78.29310 1471619 2022-06-29 NA NA
## 453  453 908374.2 76.36148 1471624 2022-06-30 NA NA

# date of peak active cases
# max_I = which(round(fitted_projected$I) == round(max(fitted_projected$I)))  # at times this works better
max_I = which(fitted_projected$I == max(fitted_projected$I))
max_date = fitted_projected$Date[max_I]
# add cumulative infected cases
fitted_projected$total_infected = fitted_projected$I + fitted_projected$R
# predicted new cases today
new_today = (fitted_projected[fitted_projected$Date == today(), ] - fitted_projected[fitted_projected$Date == today()-1, ])$total_infected
# maximum cumulative cases, date. May add to plot.
fitted_projected$Date[min(which(round(fitted_projected$total_infected) == max(round(fitted_projected$total_infected))))]

## [1] "2022-06-30"

fitted_projected[min(which(round(fitted_projected$total_infected) == max(round(fitted_projected$total_infected)))),]

##     time        S        I       R       Date  A Rm total_infected
## 453  453 908374.2 76.36148 1471624 2022-06-30 NA NA        1471701

# save for later view
write.csv(fitted_projected, "fitted_sir.csv")

# color settings
colors = c("Susceptible" = "black", "Recovered" = "green", "Infectious" = "red", "Observed Active" = "orange", "Observed Recovered" = "blue")

# plot projection data
sirplot1 = ggplot(fitted_projected, aes(x = Date)) + 
  geom_line(aes(y = I, color = "Infectious")) + 
  geom_line(aes(y = S, color = "Susceptible")) + 
  geom_line(aes(y = R, color = "Recovered")) +
  geom_point(aes(y = A, color = "Observed Active")) +
  geom_point(aes(y = Rm, color = "Observed Recovered")) +
  scale_y_continuous(labels = scales::comma) +
  scale_x_date(date_breaks = "14 day", date_labels = "%d/%m/%y") + 
  scale_colour_manual(values = colors) +
  labs(y = "Number of cases", title = paste("COVID-19 SIR model Malaysia,", name), 
       subtitle = paste("Projection from data:", start_date, "to", fitted_projected$Date[max(Day)]),
       color = paste0("Model fit:\n",
                      "R square (I) = ", round(R2_1,3), "\n",
                      "R square (R) = ", round(R2_2,3), "\n",
                      "MAPE (I) = ", round(mape1,3), "\n",
                      "MAPE (R) = ", round(mape2,3), "\n\n",
                      "SIR parameters:\n",
                      "R0 = ", round(R0, 3), "\n",
                      "beta = ", round(Opt_par[1], 3), "\n",
                      "gamma = ", round(Opt_par[2], 3), "\n\n",
                      "Susceptible = ", round(n), "\n",
                      "Peak Active = ", round(max(fitted_projected$I)), "\n",
                      "Maximum Total Infected = ", round(max(fitted_projected$total_infected)))) +
  geom_vline(xintercept = as.numeric(as.Date(max_date)), linetype = "dotted") +
  annotate(geom = "text", x = as.Date(max_date)+20, y = n*1.3, 
           label = paste0("Peak on ", format(max_date, "%d/%m/%y")), angle = 0) +
  # geom_vline(xintercept = as.numeric(as.Date(today())), linetype = "dotted", color = "red") +
  # annotate(geom = "text", x = as.Date(today())+25, y = n*1.2, 
  #          label = paste0("Today's Prediction (", format(today(), "%d/%m/%y"), ")\n",
  #                         "Total Cases = ", round(fitted_projected[fitted_projected$Date == today(), "total_infected"]),
  #                         "\nNew Cases = ", round(new_today)), angle = 0) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
sirplot1

## Warning: Removed 336 rows containing missing values (geom_point).

## Warning: Removed 336 rows containing missing values (geom_point).

# uncomment ggsave line to save plot in png format, make sure folder "plots_srn" is created
# ggsave(paste0("plots_srn/sir", name, ".png"), width = 12, height = 9)

# plot projection data, in log10
sirplot1_log = ggplot(fitted_projected, aes(x = Date)) + 
  geom_line(aes(y = I, color = "Infectious")) + 
  geom_line(aes(y = S, color = "Susceptible")) + 
  geom_line(aes(y = R, color = "Recovered")) +
  geom_point(aes(y = A, color = "Observed Active")) +
  geom_point(aes(y = Rm, color = "Observed Recovered")) +
  scale_y_log10(labels = scales::comma) +
  scale_x_date(date_breaks = "14 day", date_labels = "%d/%m/%y") + 
  scale_colour_manual(values = colors) +
  labs(y = "Number of cases", title = paste("COVID-19 SIR model Malaysia,", name, "log10"), 
       subtitle = paste("Projection from data:", start_date, "to", fitted_projected$Date[max(Day)]),
       color = paste0("Model fit:\n",
                      "R square (I) = ", round(R2_1,3), "\n",
                      "R square (R) = ", round(R2_2,3), "\n",
                      "MAPE (I) = ", round(mape1,3), "\n",
                      "MAPE (R) = ", round(mape2,3), "\n\n",
                      "SIR parameters:\n",
                      "R0 = ", round(R0, 3), "\n",
                      "beta = ", round(Opt_par[1], 3), "\n",
                      "gamma = ", round(Opt_par[2], 3), "\n\n",
                      "Susceptible = ", round(n), "\n",
                      "Peak Active = ", round(max(fitted_projected$I)), "\n",
                      "Maximum Total Infected = ", round(max(fitted_projected$total_infected)))) +
  geom_vline(xintercept = as.numeric(as.Date(max_date)), linetype = "dotted") +
  annotate(geom = "text", x = as.Date(max_date)+20, y = n*1.3, 
           label = paste0("Peak on ", format(max_date, "%d/%m/%y")), angle = 0) +
  # geom_vline(xintercept = as.numeric(as.Date(today())), linetype = "dotted", color = "red") +
  # annotate(geom = "text", x = as.Date(today())+25, y = n*0.7, 
  #          label = paste0("Today's Prediction (", format(today(), "%d/%m/%y"), ")\n",
  #                         "Total Cases = ", round(fitted_projected[fitted_projected$Date == today(), "total_infected"]),
  #                         "\nNew Cases = ", round(new_today)), angle = 0) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
sirplot1_log

## Warning: Removed 336 rows containing missing values (geom_point).

## Warning: Removed 336 rows containing missing values (geom_point).

# uncomment ggsave line to save plot in png format, make sure folder "plots_srn" is created
# ggsave(paste0("plots_srn/sir", name, "_log.png"), width = 12, height = 9)

Projection + Asymptomatic

# if 5%, 50%  asymptomatic
# https://www.cebm.net/covid-19/covid-19-what-proportion-are-asymptomatic/

# 1st situation --- #
# get the fitted values from our SIR model
p_asym = .05  # China figure, 5%
n1 = n/(1-p_asym)  # inflate by %
# init = c(S=n1-Infected[1]-Recovered[1]-Death[1], I=Infected[1]-Death[1]-Recovered[1], R=Recovered[1]+Death[1])
init_pre = init  # save original value for use later
# new init
init = init_pre/(1-p_asym)
fitted_projected1 = data.frame(ode(y=init, times=t, func=SIR, parms=Opt_par))
# add a Date & Active
fitted_projected1$Date = ymd(start_date) + days(t - 1)
fitted_projected1$A = c(Active, rep(NA, length(t) - length(Active)))
fitted_projected1$Rm = c(Removed, rep(NA, length(t) - length(Active)))
head(fitted_projected1, 10); tail(fitted_projected1, 10)

##    time       S        I        R       Date     A     Rm
## 1     1 2135912 15272.63 354157.9 2021-04-04 14509 336450
## 2     2 2134266 15829.54 355246.4 2021-04-05 14278 337751
## 3     3 2132562 16405.38 356374.5 2021-04-06 14161 339168
## 4     4 2130798 17000.71 357543.6 2021-04-07 14097 340371
## 5     5 2128971 17616.07 358755.1 2021-04-08 14203 341550
## 6     6 2127080 18252.03 360010.3 2021-04-09 14805 342802
## 7     7 2125122 18909.14 361310.8 2021-04-10 15059 344058
## 8     8 2123096 19587.98 362658.1 2021-04-11 15574 345282
## 9     9 2120999 20289.12 364053.7 2021-04-12 15835 346338
## 10   10 2118830 21013.13 365499.1 2021-04-13 16300 347640

##     time        S        I       R       Date  A Rm
## 444  444 850437.8 36.57164 1654868 2022-06-21 NA NA
## 445  445 850436.3 35.56708 1654870 2022-06-22 NA NA
## 446  446 850434.9 34.59011 1654873 2022-06-23 NA NA
## 447  447 850433.4 33.63997 1654875 2022-06-24 NA NA
## 448  448 850432.0 32.71593 1654877 2022-06-25 NA NA
## 449  449 850430.7 31.81728 1654880 2022-06-26 NA NA
## 450  450 850429.3 30.94330 1654882 2022-06-27 NA NA
## 451  451 850428.0 30.09333 1654884 2022-06-28 NA NA
## 452  452 850426.8 29.26670 1654886 2022-06-29 NA NA
## 453  453 850425.6 28.46278 1654888 2022-06-30 NA NA

# add cumulative cases
fitted_projected1$total_infected = fitted_projected1$I + fitted_projected1$R
# maximum cumulative cases, date. May add to plot.
fitted_projected1$Date[min(which(round(fitted_projected1$total_infected) == max(round(fitted_projected1$total_infected))))]

## [1] "2022-06-30"

fitted_projected1[min(which(round(fitted_projected1$total_infected) == max(round(fitted_projected1$total_infected)))),]

##     time        S        I       R       Date  A Rm total_infected
## 453  453 850425.6 28.46278 1654888 2022-06-30 NA NA        1654917

# --- 2nd situation --- #
# get the fitted values from our SIR model
p_asym = .5  # Iceland figure, 50%
n2 = n/(1-p_asym)  # inflate by %
# init = c(S=n2-Infected[1]-Recovered[1]-Death[1], I=Infected[1]-Death[1]-Recovered[1], R=Recovered[1]+Death[1])
# new init  
init = init_pre/(1-p_asym)
fitted_projected2 = data.frame(ode(y=init, times=t, func=SIR, parms=Opt_par))
# add a Date & Active
fitted_projected2$Date = ymd(start_date) + days(t - 1)
fitted_projected2$A = c(Active, rep(NA, length(t) - length(Active)))
fitted_projected2$Rm = c(Removed, rep(NA, length(t) - length(Active)))
head(fitted_projected2, 10); tail(fitted_projected2, 10)

##    time       S        I        R       Date     A     Rm
## 1     1 4058232 29018.00 672900.0 2021-04-04 14509 336450
## 2     2 4052001 33078.71 675070.3 2021-04-05 14278 337751
## 3     3 4044911 37695.27 677544.0 2021-04-06 14161 339168
## 4     4 4036848 42940.00 680362.3 2021-04-07 14097 340371
## 5     5 4027684 48893.59 683572.1 2021-04-08 14203 341550
## 6     6 4017278 55645.71 687226.1 2021-04-09 14805 342802
## 7     7 4005471 63295.47 691383.5 2021-04-10 15059 344058
## 8     8 3992087 71951.99 696110.9 2021-04-11 15574 345282
## 9     9 3976932 81734.64 701483.0 2021-04-12 15835 346338
## 10   10 3959794 92773.17 707583.1 2021-04-13 16300 347640

##     time        S            I       R       Date  A Rm
## 444  444 272512.1 0.0003857938 4487638 2022-06-21 NA NA
## 445  445 272512.1 0.0003646029 4487638 2022-06-22 NA NA
## 446  446 272512.1 0.0003445759 4487638 2022-06-23 NA NA
## 447  447 272512.1 0.0003256489 4487638 2022-06-24 NA NA
## 448  448 272512.1 0.0003077616 4487638 2022-06-25 NA NA
## 449  449 272512.1 0.0002908568 4487638 2022-06-26 NA NA
## 450  450 272512.1 0.0002748806 4487638 2022-06-27 NA NA
## 451  451 272512.1 0.0002597819 4487638 2022-06-28 NA NA
## 452  452 272512.1 0.0002455125 4487638 2022-06-29 NA NA
## 453  453 272512.1 0.0002320270 4487638 2022-06-30 NA NA

# add cumulative cases
fitted_projected2$total_infected = fitted_projected2$I + fitted_projected2$R
# maximum cumulative cases, date. May add to plot.
fitted_projected2$Date[min(which(round(fitted_projected2$total_infected) == max(round(fitted_projected2$total_infected))))]

## [1] "2022-01-25"

fitted_projected2[min(which(round(fitted_projected2$total_infected) == max(round(fitted_projected2$total_infected)))),]

##     time        S        I       R       Date  A Rm total_infected
## 297  297 272512.5 1.559656 4487636 2022-01-25 NA NA        4487638

# color settings
colors = c("Susceptible" = "black", "Recovered" = "green", "Infectious" = "red", 
           "Susceptible + 5% (dashed)" = "black", "Recovered + 5% (dashed)" = "green", "Infectious + 5% (dashed)" = "red", 
           "Susceptible + 50% (dot-dash)" = "black", "Recovered + 50% (dot-dash)" = "green", "Infectious + 50% (dot-dash)" = "red", 
           "Observed Active" = "orange", "Observed Recovered" = "blue")

# plot projection data
sirplot2 = ggplot(fitted_projected, aes(x = Date)) + 
  geom_line(aes(y = I, color = "Infectious")) + 
  geom_line(aes(y = S, color = "Susceptible")) + 
  geom_line(aes(y = R, color = "Recovered")) +
  geom_line(data=fitted_projected1, aes(y = I, color = "Infectious + 5% (dashed)"), linetype = "dashed") + 
  geom_line(data=fitted_projected1, aes(y = S, color = "Susceptible + 5% (dashed)"), linetype = "dashed") + 
  geom_line(data=fitted_projected1, aes(y = R, color = "Recovered + 5% (dashed)"), linetype = "dashed") +
  geom_line(data=fitted_projected2, aes(y = I, color = "Infectious + 50% (dot-dash)"), linetype = "dotdash") + 
  geom_line(data=fitted_projected2, aes(y = S, color = "Susceptible + 50% (dot-dash)"), linetype = "dotdash") + 
  geom_line(data=fitted_projected2, aes(y = R, color = "Recovered + 50% (dot-dash)"), linetype = "dotdash") +
  geom_point(aes(y = A, color = "Observed Active")) +
  geom_point(aes(y = Rm, color = "Observed Recovered")) +
  scale_y_continuous(labels = scales::comma) +
  scale_x_date(date_breaks = "14 day", date_labels = "%d/%m/%y") + 
  scale_colour_manual(values = colors) +
  labs(y = "Number of cases", title = paste("COVID-19 SIR model Malaysia,", name), 
       subtitle = paste("Projection from data:", start_date, "to", fitted_projected$Date[max(Day)]),
       color = paste0("Model fit:\n",
                      "R square (I) = ", round(R2_1,3), "\n",
                      "R square (R) = ", round(R2_2,3), "\n",
                      "MAPE (I) = ", round(mape1,3), "\n",
                      "MAPE (R) = ", round(mape2,3), "\n\n",
                      "SIR parameters:\n",
                      "R0 = ", round(R0, 3), "\n",
                      "beta = ", round(Opt_par[1], 3), "\n",
                      "gamma = ", round(Opt_par[2], 3), "\n\n",
                      "Susceptible Individuals:\n",
                      "No asymptomatic = ", round(n), "\n",
                      "Asymptomatic 5% = ", round(n1), "\n",
                      "Asymptomatic 50% = ", round(n2), "\n\n",
                      "Peak Active Cases:\n",
                      "No asymptomatic = ", round(max(fitted_projected$I)), "\n",
                      "Asymptomatic 5% = ", round(max(fitted_projected1$I)), "\n",
                      "Asymptomatic 50% = ", round(max(fitted_projected2$I)), "\n\n",
                      "Maximum Total Infected:\n",
                      "No asymptomatic = ", round(max(fitted_projected$total_infected)), "\n",
                      "Asymptomatic 5% = ", round(max(fitted_projected1$total_infected)), "\n",
                      "Asymptomatic 50% = ", round(max(fitted_projected2$total_infected)), "\n")) +
  geom_vline(xintercept = as.numeric(as.Date(max_date)), linetype = "dotted") +
  annotate(geom = "text", x = as.Date(max_date)+20, y = n2*1.3, 
           label = paste0("Peak on ", format(max_date, "%d/%m/%y")), angle = 0) +
  # geom_vline(xintercept = as.numeric(as.Date(today())), linetype = "dotted", color = "red") +
  # annotate(geom = "text", x = as.Date(today())+25, y = n2*1.2, 
  #          label = paste0("Today's Prediction (", format(today(), "%d/%m/%y"), ")\n",
  #                         "Total Cases = ", round(fitted_projected[fitted_projected$Date == today(), "total_infected"]),
  #                         "\nNew Cases = ", round(new_today)), angle = 0) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
sirplot2

## Warning: Removed 336 rows containing missing values (geom_point).

## Warning: Removed 336 rows containing missing values (geom_point).

# uncomment ggsave line to save plot in png format, make sure folder "plots_srn" is created
ggsave(paste0("plots_srn/sir", name, "_", fitted_projected$Date[max(Day)], "_asymp.png"), width = 12, height = 9)

## Warning: Removed 336 rows containing missing values (geom_point).

## Warning: Removed 336 rows containing missing values (geom_point).

# plot projection data, in log10
sirplot2_log = ggplot(fitted_projected, aes(x = Date)) + 
  geom_line(aes(y = I, color = "Infectious")) + 
  geom_line(aes(y = S, color = "Susceptible")) + 
  geom_line(aes(y = R, color = "Recovered")) +
  geom_line(data=fitted_projected1, aes(y = I, color = "Infectious + 5% (dashed)"), linetype = "dashed") + 
  geom_line(data=fitted_projected1, aes(y = S, color = "Susceptible + 5% (dashed)"), linetype = "dashed") + 
  geom_line(data=fitted_projected1, aes(y = R, color = "Recovered + 5% (dashed)"), linetype = "dashed") +
  geom_line(data=fitted_projected2, aes(y = I, color = "Infectious + 50% (dot-dash)"), linetype = "dotdash") + 
  geom_line(data=fitted_projected2, aes(y = S, color = "Susceptible + 50% (dot-dash)"), linetype = "dotdash") + 
  geom_line(data=fitted_projected2, aes(y = R, color = "Recovered + 50% (dot-dash)"), linetype = "dotdash") +
  geom_point(aes(y = A, color = "Observed Active")) +
  geom_point(aes(y = Rm, color = "Observed Recovered")) +
  scale_y_log10(labels = scales::comma) +
  scale_x_date(date_breaks = "14 day", date_labels = "%d/%m/%y") + 
  scale_colour_manual(values = colors) +
  labs(y = "Number of cases", title = paste("COVID-19 SIR model Malaysia,", name), 
       subtitle = paste("Projection from data:", start_date, "to", fitted_projected$Date[max(Day)]),
       color = paste0("Model fit:\n",
                      "R square (I) = ", round(R2_1,3), "\n",
                      "R square (R) = ", round(R2_2,3), "\n",
                      "MAPE (I) = ", round(mape1,3), "\n",
                      "MAPE (R) = ", round(mape2,3), "\n\n",
                      "SIR parameters:\n",
                      "R0 = ", round(R0, 3), "\n",
                      "beta = ", round(Opt_par[1], 3), "\n",
                      "gamma = ", round(Opt_par[2], 3), "\n\n",
                      "Susceptible Individuals:\n",
                      "No asymptomatic = ", round(n), "\n",
                      "Asymptomatic 5% = ", round(n1), "\n",
                      "Asymptomatic 50% = ", round(n2), "\n\n",
                      "Peak Active Cases:\n",
                      "No asymptomatic = ", round(max(fitted_projected$I)), "\n",
                      "Asymptomatic 5% = ", round(max(fitted_projected1$I)), "\n",
                      "Asymptomatic 50% = ", round(max(fitted_projected2$I)), "\n\n",
                      "Maximum Total Infected:\n",
                      "No asymptomatic = ", round(max(fitted_projected$total_infected)), "\n",
                      "Asymptomatic 5% = ", round(max(fitted_projected1$total_infected)), "\n",
                      "Asymptomatic 50% = ", round(max(fitted_projected2$total_infected)), "\n")) +
  geom_vline(xintercept = as.numeric(as.Date(max_date)), linetype = "dotted") +
  annotate(geom = "text", x = as.Date(max_date)+20, y = n2*1.3, 
           label = paste0("Peak on ", format(max_date, "%d/%m/%y")), angle = 0) +
  # geom_vline(xintercept = as.numeric(as.Date(today())), linetype = "dotted", color = "red") +
  # annotate(geom = "text", x = as.Date(today())+25, y = n2*0.7, 
  #          label = paste0("Today's Prediction (", format(today(), "%d/%m/%y"), ")\n",
  #                         "Total Cases = ", round(fitted_projected[fitted_projected$Date == today(), "total_infected"]),
  #                         "\nNew Cases = ", round(new_today)), angle = 0) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
sirplot2_log

## Warning: Removed 336 rows containing missing values (geom_point).

## Warning: Removed 336 rows containing missing values (geom_point).

# uncomment ggsave line to save plot in png format, make sure folder "plots_srn" is created
ggsave(paste0("plots_srn/sir", name, "_", fitted_projected$Date[max(Day)], "_log_asymp.png"), width = 12, height = 9)

## Warning: Removed 336 rows containing missing values (geom_point).

## Warning: Removed 336 rows containing missing values (geom_point).

Appendix

Detailed projected data.

fitted_projected$new_cases = diff(c(NA, fitted_projected$total_infected))
knitr::kable(fitted_projected, format="markdown")

time	S	I	R	Date	A	Rm	total_infected	new_cases
1	2029116.0	14509.00000	336450.0	2021-04-04	14509	336450	350959.0	NA
2	2027634.9	14958.73911	337481.3	2021-04-05	14278	337751	352440.1	1481.064297
3	2026109.2	15421.26667	338544.6	2021-04-06	14161	339168	353965.8	1525.771812
4	2024537.5	15896.86238	339640.6	2021-04-07	14097	340371	355537.5	1571.633247
5	2022918.7	16385.83710	340770.4	2021-04-08	14203	341550	357156.3	1618.791162
6	2021251.6	16888.47719	341934.9	2021-04-09	14805	342802	358823.4	1667.161773
7	2019534.8	17405.07638	343135.1	2021-04-10	15059	344058	360540.2	1716.791678
8	2017767.1	17935.92907	344372.0	2021-04-11	15574	345282	362307.9	1767.703458
9	2015947.2	18481.32973	345646.5	2021-04-12	15835	346338	364127.8	1819.918865
10	2014073.7	19041.57241	346959.7	2021-04-13	16300	347640	366001.3	1873.456816
11	2012145.4	19616.95023	348312.7	2021-04-14	16696	349133	367929.6	1928.336963
12	2010160.8	20207.75483	349706.5	2021-04-15	17575	350402	369914.2	1984.578442
13	2008118.6	20814.27585	351142.1	2021-04-16	18600	351928	371956.4	2042.199575
14	2006017.4	21436.80018	352620.8	2021-04-17	19094	353765	374057.6	2101.217811
15	2003855.7	22075.61138	354143.7	2021-04-18	19854	355200	376219.3	2161.649777
16	2001632.2	22730.98896	355711.8	2021-04-19	20522	356610	378442.8	2223.511145
17	1999345.4	23403.20762	357326.4	2021-04-20	21268	358205	380729.6	2286.816510
18	1996993.8	24092.53648	358988.6	2021-04-21	21687	360126	383081.2	2351.579318
19	1994576.0	24799.23823	360699.8	2021-04-22	22014	362674	385499.0	2417.811772
20	1992090.5	25523.56829	362461.0	2021-04-23	22512	365023	387984.5	2485.524718
21	1989535.8	26265.77384	364273.5	2021-04-24	22926	367326	390539.2	2554.727549
22	1986910.3	27026.09293	366138.6	2021-04-25	23753	369189	393164.7	2625.428092
23	1984212.7	27804.75344	368057.6	2021-04-26	24713	371005	395862.3	2697.632498
24	1981441.3	28601.97204	370031.7	2021-04-27	25414	373037	398633.7	2771.345126
25	1978594.8	29417.95311	372062.3	2021-04-28	26719	374874	401480.2	2846.568424
26	1975671.5	30252.88759	374150.6	2021-04-29	28093	376832	404403.5	2923.302812
27	1972669.9	31106.95184	376298.1	2021-04-30	29227	379486	407405.1	3001.546558
28	1969588.6	31980.30642	378506.1	2021-05-01	29631	381963	410486.4	3081.295655
29	1966426.1	32873.09482	380775.8	2021-05-02	30339	384673	413648.9	3162.543695
30	1963180.8	33785.44219	383108.7	2021-05-03	30753	386759	416894.2	3245.281743
31	1959851.3	34717.45401	385506.2	2021-05-04	31516	389116	420223.7	3329.498209
32	1956436.1	35669.21474	387969.7	2021-05-05	32939	391437	423638.9	3415.178722
33	1952933.8	36640.78640	390500.4	2021-05-06	33762	394165	427141.2	3502.306005
34	1949343.0	37632.20722	393099.8	2021-05-07	34789	397636	430732.0	3590.859743
35	1945662.1	38643.49013	395769.4	2021-05-08	36564	400380	434412.9	3680.816468
36	1941890.0	39674.62135	398510.4	2021-05-09	37060	403617	438185.0	3772.149430
37	1938025.2	40725.55889	401324.3	2021-05-10	37396	407088	442049.8	3864.828483
38	1934066.4	41796.23106	404212.4	2021-05-11	38499	409958	446008.6	3958.819969
39	1930012.3	42886.53501	407176.2	2021-05-12	40101	413121	450062.7	4054.086606
40	1925861.7	43996.33521	410217.0	2021-05-13	41582	416495	454213.3	4150.587389
41	1921613.4	45125.46198	413336.1	2021-05-14	41471	420719	458461.6	4248.277485
42	1917266.3	46273.70998	416535.0	2021-05-15	42135	424195	462808.7	4347.108149
43	1912819.3	47440.83683	419814.9	2021-05-16	41889	428221	467255.7	4447.026639
44	1908271.3	48626.56162	423177.1	2021-05-17	43506	431050	471803.7	4547.976137
45	1903621.4	49830.56355	426623.0	2021-05-18	44827	434594	476453.6	4649.895698
46	1898868.7	51052.48056	430153.8	2021-05-19	47340	438156	481206.3	4752.720185
47	1894012.3	52291.90808	433770.8	2021-05-20	50171	442131	486062.7	4856.380240
48	1889051.5	53548.39777	437475.1	2021-05-21	52106	446689	491023.5	4960.802252
49	1883985.6	54821.45635	441268.0	2021-05-22	53682	451433	496089.4	5065.908346
50	1878814.0	56110.54456	445150.5	2021-05-23	57022	455069	501261.0	5171.616386
51	1873536.1	57415.07612	449123.8	2021-05-24	60018	458582	506538.9	5277.839995
52	1868151.6	58734.41690	453188.9	2021-05-25	63458	462431	511923.4	5384.488588
53	1862660.2	60067.88412	457346.9	2021-05-26	66208	467159	517414.8	5491.467431
54	1857061.5	61414.74569	461598.8	2021-05-27	69408	471816	523013.5	5598.677709
55	1851355.5	62774.21972	465945.3	2021-05-28	72823	476691	528719.5	5706.016620
56	1845542.1	64145.47412	470387.4	2021-05-29	76218	482316	534532.9	5813.377487
57	1839621.4	65527.62640	474925.9	2021-05-30	78017	487516	540453.6	5920.649890
58	1833593.7	66919.74358	479561.5	2021-05-31	79523	492834	546481.3	6027.719819
59	1827459.3	68320.84236	484294.9	2021-06-01	80474	498988	552615.7	6134.469854
60	1821218.5	69729.88936	489126.6	2021-06-02	82274	504891	558856.5	6240.779356
61	1814872.0	71145.80166	494057.2	2021-06-03	83331	512043	565203.0	6346.524683
62	1808420.4	72567.44743	499087.2	2021-06-04	84369	518753	571654.6	6451.579437
63	1801864.6	73993.64689	504216.8	2021-06-05	85607	524967	578210.4	6555.814718
64	1795205.5	75423.17333	509446.4	2021-06-06	86628	530187	584869.5	6659.099412
65	1788444.2	76854.75454	514776.1	2021-06-07	84269	537817	591630.8	6761.300486
66	1781581.9	78287.07424	520206.0	2021-06-08	82797	544855	598493.1	6862.283315
67	1774620.0	79718.77392	525736.3	2021-06-09	81575	552316	605455.0	6961.912020
68	1767559.9	81148.45481	531366.6	2021-06-10	79848	559714	612515.1	7060.049823
69	1760403.4	82574.68009	537097.0	2021-06-11	78864	567547	619671.6	7156.559426
70	1753152.1	83995.97737	542927.0	2021-06-12	76247	575957	626922.9	7251.303393
71	1745807.9	85410.84130	548856.3	2021-06-13	73324	584184	634267.1	7344.144557
72	1738373.0	86817.73649	554884.3	2021-06-14	71625	590832	641702.0	7434.946428
73	1730849.4	88215.10061	561010.5	2021-06-15	70112	597764	649225.6	7523.573619
74	1723239.5	89601.34768	567234.2	2021-06-16	67949	605077	656835.5	7609.892272
75	1715545.7	90974.87156	573554.4	2021-06-17	66097	612667	664529.3	7693.770494
76	1707770.6	92334.04967	579970.3	2021-06-18	65602	619602	672304.4	7775.078792
77	1699917.0	93677.24680	586480.8	2021-06-19	64523	626592	680158.0	7853.690513
78	1691987.5	95002.81915	593084.7	2021-06-20	63815	632593	688087.5	7929.482276
79	1683985.1	96309.11850	599780.7	2021-06-21	62918	638101	696089.9	8002.334403
80	1675913.0	97594.49647	606567.5	2021-06-22	62027	643735	704162.0	8072.131341
81	1667774.2	98857.30895	613443.4	2021-06-23	60816	650190	712300.8	8138.762076
82	1659572.1	100095.92052	620407.0	2021-06-24	61162	655685	720502.9	8202.120530
83	1651310.0	101308.70910	627456.3	2021-06-25	60117	662542	728765.0	8262.105947
84	1642991.4	102494.07045	634589.5	2021-06-26	60646	667816	737083.6	8318.623262
85	1634619.8	103650.42284	641804.8	2021-06-27	61395	672653	745455.2	8371.583444
86	1626198.9	104776.21169	649099.9	2021-06-28	61812	677454	753876.1	8420.903822
87	1617732.4	105869.91412	656472.7	2021-06-29	62844	682859	762342.6	8466.508384
88	1609224.1	106930.04355	663920.9	2021-06-30	64129	687850	770850.9	8508.328053
89	1600677.8	107955.15415	671442.1	2021-07-01	65453	693514	779397.2	8546.300930
90	1592097.4	108943.84524	679033.8	2021-07-02	66084	699865	787977.6	8580.372509
91	1583486.9	109894.76554	686693.3	2021-07-03	66958	705649	796588.1	8610.495859
92	1574850.3	110806.61734	694418.1	2021-07-04	67669	710983	805224.7	8636.631782
93	1566191.5	111678.16041	702205.3	2021-07-05	69447	715592	813883.5	8658.748923
94	1557514.7	112508.21577	710052.1	2021-07-06	72201	720492	822560.3	8676.823859
95	1548823.9	113295.66930	717955.5	2021-07-07	74344	725446	831251.1	8690.841147
96	1540123.1	114039.47505	725912.5	2021-07-08	77275	731383	839951.9	8700.793340
97	1531416.4	114738.65835	733920.0	2021-07-09	80665	737173	848658.6	8706.680965
98	1522707.9	115392.31859	741974.8	2021-07-10	84021	743170	857367.1	8708.512474
99	1514001.6	115999.63186	750073.8	2021-07-11	87841	748455	866073.4	8706.304152
100	1505301.5	116559.85316	758213.7	2021-07-12	91272	753598	874773.5	8700.080003
101	1496611.6	117072.31839	766391.1	2021-07-13	96236	759713	883463.4	8689.871592
102	1487935.9	117536.44599	774602.7	2021-07-14	101359	766208	892139.1	8675.717874
103	1479278.2	117951.73835	782845.0	2021-07-15	108369	772413	900796.8	8657.664973
104	1470642.5	118317.78274	791114.8	2021-07-16	114053	779270	909432.5	8635.765957
105	1462032.4	118634.25210	799408.4	2021-07-17	119814	786037	918042.6	8610.080566
106	1453451.7	118900.90539	807722.4	2021-07-18	124593	791968	926623.3	8580.674934
107	1444904.1	119117.58763	816053.3	2021-07-19	128997	798536	935170.9	8547.621284
108	1436393.1	119284.22968	824397.7	2021-07-20	133703	806196	943681.9	8510.997600
109	1427922.2	119400.84765	832751.9	2021-07-21	137587	814297	952152.8	8470.887288
110	1419494.8	119467.54207	841112.6	2021-07-22	142051	822867	960580.2	8427.378824
111	1411114.3	119484.49664	849476.2	2021-07-23	147386	833105	968960.7	8380.565386
112	1402783.7	119451.97686	857839.3	2021-07-24	153633	842760	977291.3	8330.544477
113	1394506.3	119370.32830	866198.4	2021-07-25	160903	852535	985568.7	8277.417546
114	1386285.0	119239.97456	874550.0	2021-07-26	165840	862114	993790.0	8221.289601
115	1378122.7	119061.41513	882890.8	2021-07-27	170224	873847	1001952.3	8162.268815
116	1370022.3	118835.22289	891217.5	2021-07-28	175113	886363	1010052.7	8100.466144
117	1361986.3	118562.04147	899526.7	2021-07-29	179179	899467	1018088.7	8035.994935
118	1354017.3	118242.58242	907815.1	2021-07-30	NA	NA	1026057.7	7968.970541
119	1346117.8	117877.62217	916079.6	2021-07-31	NA	NA	1033957.2	7899.509951
120	1338290.1	117467.99886	924316.9	2021-08-01	NA	NA	1041784.9	7827.731410
121	1330536.3	117014.60906	932524.1	2021-08-02	NA	NA	1049538.7	7753.754066
122	1322858.6	116518.40433	940698.0	2021-08-03	NA	NA	1057216.4	7677.697616
123	1315258.9	115980.38769	948835.7	2021-08-04	NA	NA	1064816.1	7599.681969
124	1307739.1	115401.61007	956934.3	2021-08-05	NA	NA	1072335.9	7519.826922
125	1300300.9	114783.16662	964991.0	2021-08-06	NA	NA	1079774.1	7438.251846
126	1292945.8	114126.19305	973003.0	2021-08-07	NA	NA	1087129.2	7355.075394
127	1285675.4	113431.86190	980967.8	2021-08-08	NA	NA	1094399.6	7270.415217
128	1278491.0	112701.37887	988882.6	2021-08-09	NA	NA	1101584.0	7184.387705
129	1271393.9	111935.97905	996745.2	2021-08-10	NA	NA	1108681.1	7097.107735
130	1264385.2	111136.92330	1004552.9	2021-08-11	NA	NA	1115689.8	7008.688445
131	1257465.9	110305.49463	1012303.6	2021-08-12	NA	NA	1122609.1	6919.241022
132	1250637.1	109442.99459	1019994.9	2021-08-13	NA	NA	1129437.9	6828.874503
133	1243899.4	108550.73978	1027624.9	2021-08-14	NA	NA	1136175.6	6737.695601
134	1237253.6	107630.05846	1035191.4	2021-08-15	NA	NA	1142821.4	6645.808541
135	1230700.2	106682.28716	1042692.5	2021-08-16	NA	NA	1149374.8	6553.314920
136	1224239.9	105708.76755	1050126.3	2021-08-17	NA	NA	1155835.1	6460.313578
137	1217873.0	104710.84321	1057491.1	2021-08-18	NA	NA	1162202.0	6366.900485
138	1211599.9	103689.85674	1064785.3	2021-08-19	NA	NA	1168475.1	6273.168649
139	1205420.7	102647.14683	1072007.2	2021-08-20	NA	NA	1174654.3	6179.208034
140	1199335.5	101584.04551	1079155.4	2021-08-21	NA	NA	1180739.5	6085.105494
141	1193344.6	100501.87560	1086228.5	2021-08-22	NA	NA	1186730.4	5990.944721
142	1187447.8	99401.94818	1093225.3	2021-08-23	NA	NA	1192627.2	5896.806208
143	1181645.0	98285.56034	1100144.4	2021-08-24	NA	NA	1198430.0	5802.767221
144	1175936.1	97153.99292	1106984.9	2021-08-25	NA	NA	1204138.9	5708.901786
145	1170320.8	96008.50852	1113745.6	2021-08-26	NA	NA	1209754.2	5615.280686
146	1164798.9	94850.34960	1120425.8	2021-08-27	NA	NA	1215276.1	5521.971465
147	1159369.8	93680.73673	1127024.4	2021-08-28	NA	NA	1220705.2	5429.038448
148	1154033.3	92500.86698	1133540.8	2021-08-29	NA	NA	1226041.7	5336.542768
149	1148788.8	91311.91247	1139974.3	2021-08-30	NA	NA	1231286.2	5244.542395
150	1143635.7	90115.01903	1146324.3	2021-08-31	NA	NA	1236439.3	5153.092180
151	1138573.4	88911.30502	1152590.3	2021-09-01	NA	NA	1241501.6	5062.243905
152	1133601.4	87701.86026	1158771.8	2021-09-02	NA	NA	1246473.6	4972.046331
153	1128718.8	86487.74513	1164868.4	2021-09-03	NA	NA	1251356.2	4882.545265
154	1123925.0	85269.98974	1170880.0	2021-09-04	NA	NA	1256150.0	4793.783619
155	1119219.2	84049.59327	1176806.2	2021-09-05	NA	NA	1260855.8	4705.801481
156	1114600.6	82827.52337	1182646.9	2021-09-06	NA	NA	1265474.4	4618.636186
157	1110068.3	81604.71576	1188402.0	2021-09-07	NA	NA	1270006.7	4532.322391
158	1105621.4	80382.07382	1194071.5	2021-09-08	NA	NA	1274453.6	4446.892157
159	1101259.0	79160.46839	1199655.5	2021-09-09	NA	NA	1278816.0	4362.375022
160	1096980.2	77940.73759	1205154.0	2021-09-10	NA	NA	1283094.8	4278.798088
161	1092784.0	76723.68676	1210567.3	2021-09-11	NA	NA	1287291.0	4196.186102
162	1088669.5	75510.08850	1215895.4	2021-09-12	NA	NA	1291405.5	4114.561541
163	1084635.5	74300.68276	1221138.8	2021-09-13	NA	NA	1295439.5	4033.944695
164	1080681.2	73096.17700	1226297.7	2021-09-14	NA	NA	1299393.8	3954.353750
165	1076805.4	71897.24649	1231372.4	2021-09-15	NA	NA	1303269.6	3875.804876
166	1073007.1	70704.53455	1236363.4	2021-09-16	NA	NA	1307067.9	3798.312308
167	1069285.2	69518.65299	1241271.2	2021-09-17	NA	NA	1310789.8	3721.888429
168	1065638.6	68340.18247	1246096.2	2021-09-18	NA	NA	1314436.4	3646.543854
169	1062066.3	67169.67303	1250839.0	2021-09-19	NA	NA	1318008.7	3572.287509
170	1058567.2	66007.64456	1255500.1	2021-09-20	NA	NA	1321507.8	3499.126713
171	1055140.1	64854.58740	1260080.3	2021-09-21	NA	NA	1324934.9	3427.067252
172	1051784.0	63710.96292	1264580.0	2021-09-22	NA	NA	1328291.0	3356.113464
173	1048497.8	62577.20414	1269000.0	2021-09-23	NA	NA	1331577.2	3286.268303
174	1045280.2	61453.71638	1273341.1	2021-09-24	NA	NA	1334794.8	3217.533420
175	1042130.3	60340.87799	1277603.8	2021-09-25	NA	NA	1337944.7	3149.909234
176	1039046.9	59239.04098	1281789.0	2021-09-26	NA	NA	1341028.1	3083.394996
177	1036028.9	58148.53178	1285897.5	2021-09-27	NA	NA	1344046.1	3017.988861
178	1033075.2	57069.65197	1289930.1	2021-09-28	NA	NA	1346999.8	2953.687951
179	1030184.8	56002.67902	1293887.6	2021-09-29	NA	NA	1349890.2	2890.488418
180	1027356.4	54947.86706	1297770.8	2021-09-30	NA	NA	1352718.6	2828.385507
181	1024589.0	53905.44758	1301580.6	2021-10-01	NA	NA	1355486.0	2767.373611
182	1021881.5	52875.63027	1305317.8	2021-10-02	NA	NA	1358193.5	2707.446331
183	1019233.0	51858.60373	1308983.4	2021-10-03	NA	NA	1360842.0	2648.596527
184	1016642.1	50854.53625	1312578.3	2021-10-04	NA	NA	1363432.9	2590.816373
185	1014108.0	49863.57657	1316103.4	2021-10-05	NA	NA	1365967.0	2534.097408
186	1011629.6	48885.85463	1319559.5	2021-10-06	NA	NA	1368445.4	2478.430581
187	1009205.8	47921.48230	1322947.7	2021-10-07	NA	NA	1370869.2	2423.806297
188	1006835.6	46970.55416	1326268.9	2021-10-08	NA	NA	1373239.4	2370.214465
189	1004517.9	46033.14819	1329523.9	2021-10-09	NA	NA	1375557.1	2317.644535
190	1002251.9	45109.32651	1332713.8	2021-10-10	NA	NA	1377823.1	2266.085541
191	1000036.3	44199.13608	1335839.5	2021-10-11	NA	NA	1380038.7	2215.526141
192	997870.4	43302.60942	1338902.0	2021-10-12	NA	NA	1382204.6	2165.954649
193	995753.0	42419.76522	1341902.2	2021-10-13	NA	NA	1384322.0	2117.359071
194	993683.3	41550.60910	1344841.1	2021-10-14	NA	NA	1386391.7	2069.727139
195	991660.2	40695.13419	1347719.6	2021-10-15	NA	NA	1388414.8	2023.046343
196	989682.9	39853.32181	1350538.7	2021-10-16	NA	NA	1390392.1	1977.303958
197	987750.5	39025.14206	1353299.4	2021-10-17	NA	NA	1392324.5	1932.487071
198	985861.9	38210.55445	1356002.6	2021-10-18	NA	NA	1394213.1	1888.582614
199	984016.3	37409.50846	1358649.2	2021-10-19	NA	NA	1396058.7	1845.577380
200	982212.8	36621.94412	1361240.2	2021-10-20	NA	NA	1397862.2	1803.458053
201	980450.6	35847.79259	1363776.6	2021-10-21	NA	NA	1399624.4	1762.211229
202	978728.8	35086.97665	1366259.2	2021-10-22	NA	NA	1401346.2	1721.823433
203	977046.5	34339.41125	1368689.1	2021-10-23	NA	NA	1403028.5	1682.281144
204	975402.9	33605.00400	1371067.0	2021-10-24	NA	NA	1404672.1	1643.570809
205	973797.3	32883.65565	1373394.1	2021-10-25	NA	NA	1406277.7	1605.678866
206	972228.7	32175.26058	1375671.1	2021-10-26	NA	NA	1407846.3	1568.591751
207	970696.4	31479.70720	1377898.9	2021-10-27	NA	NA	1409378.6	1532.295922
208	969199.6	30796.87845	1380078.5	2021-10-28	NA	NA	1410875.4	1496.777868
209	967737.6	30126.65215	1382210.8	2021-10-29	NA	NA	1412337.4	1462.024125
210	966309.6	29468.90145	1384296.5	2021-10-30	NA	NA	1413765.4	1428.021286
211	964914.8	28823.49520	1386336.7	2021-10-31	NA	NA	1415160.2	1394.756012
212	963552.6	28190.29830	1388332.1	2021-11-01	NA	NA	1416522.4	1362.215045
213	962222.2	27569.17205	1390283.6	2021-11-02	NA	NA	1417852.8	1330.385216
214	960922.9	26959.97454	1392192.1	2021-11-03	NA	NA	1419152.1	1299.253455
215	959654.1	26362.56091	1394058.3	2021-11-04	NA	NA	1420420.9	1268.806796
216	958415.1	25776.78368	1395883.1	2021-11-05	NA	NA	1421659.9	1239.032391
217	957205.2	25202.49307	1397667.3	2021-11-06	NA	NA	1422869.8	1209.917508
218	956023.7	24639.53723	1399411.7	2021-11-07	NA	NA	1424051.3	1181.449548
219	954870.1	24087.76255	1401117.1	2021-11-08	NA	NA	1425204.9	1153.616040
220	953743.7	23547.01389	1402784.3	2021-11-09	NA	NA	1426331.3	1126.404653
221	952643.9	23017.13483	1404413.9	2021-11-10	NA	NA	1427431.1	1099.803200
222	951570.1	22497.96790	1406006.9	2021-11-11	NA	NA	1428504.9	1073.799640
223	950521.7	21989.35479	1407563.9	2021-11-12	NA	NA	1429553.3	1048.382083
224	949498.2	21491.13656	1409085.7	2021-11-13	NA	NA	1430576.8	1023.538792
225	948498.9	21003.15383	1410572.9	2021-11-14	NA	NA	1431576.1	999.258189
226	947523.4	20525.24696	1412026.3	2021-11-15	NA	NA	1432551.6	975.528853
227	946571.1	20057.25626	1413446.7	2021-11-16	NA	NA	1433503.9	952.339528
228	945641.4	19599.02208	1414834.6	2021-11-17	NA	NA	1434433.6	929.679117
229	944733.9	19150.38504	1416190.8	2021-11-18	NA	NA	1435341.1	907.536692
230	943848.0	18711.18613	1417515.9	2021-11-19	NA	NA	1436227.0	885.901488
231	942983.2	18281.26686	1418810.5	2021-11-20	NA	NA	1437091.8	864.762910
232	942139.1	17860.46937	1420075.4	2021-11-21	NA	NA	1437935.9	844.110527
233	941315.1	17448.63658	1421311.2	2021-11-22	NA	NA	1438759.9	823.934080
234	940510.9	17045.61227	1422518.5	2021-11-23	NA	NA	1439564.1	804.223475
235	939726.0	16651.24121	1423697.8	2021-11-24	NA	NA	1440349.0	784.968789
236	938959.8	16265.36924	1424849.8	2021-11-25	NA	NA	1441115.2	766.160265
237	938212.0	15887.84336	1425975.2	2021-11-26	NA	NA	1441863.0	747.788313
238	937482.2	15518.51180	1427074.3	2021-11-27	NA	NA	1442592.8	729.843513
239	936769.8	15157.22414	1428147.9	2021-11-28	NA	NA	1443305.2	712.316609
240	936074.6	14803.83133	1429196.5	2021-11-29	NA	NA	1444000.4	695.198510
241	935396.2	14458.18579	1430220.6	2021-11-30	NA	NA	1444678.8	678.480292
242	934734.0	14120.14144	1431220.8	2021-12-01	NA	NA	1445341.0	662.153193
243	934087.8	13789.55380	1432197.6	2021-12-02	NA	NA	1445987.2	646.208610
244	933457.2	13466.27997	1433151.6	2021-12-03	NA	NA	1446617.8	630.638106
245	932841.7	13150.17875	1434083.1	2021-12-04	NA	NA	1447233.3	615.433400
246	932241.1	12841.11061	1434992.7	2021-12-05	NA	NA	1447833.9	600.586370
247	931655.1	12538.93779	1435881.0	2021-12-06	NA	NA	1448419.9	586.089048
248	931083.1	12243.52428	1436748.4	2021-12-07	NA	NA	1448991.9	571.933623
249	930525.0	11954.73587	1437595.3	2021-12-08	NA	NA	1449550.0	558.112435
250	929980.4	11672.44017	1438422.2	2021-12-09	NA	NA	1450094.6	544.617976
251	929449.0	11396.50665	1439229.5	2021-12-10	NA	NA	1450626.0	531.442886
252	928930.4	11126.80663	1440017.8	2021-12-11	NA	NA	1451144.6	518.579952
253	928424.3	10863.21330	1440787.4	2021-12-12	NA	NA	1451650.7	506.022109
254	927930.6	10605.60174	1441538.8	2021-12-13	NA	NA	1452144.4	493.762430
255	927448.8	10353.84892	1442272.4	2021-12-14	NA	NA	1452626.2	481.794135
256	926978.7	10107.83372	1442988.5	2021-12-15	NA	NA	1453096.3	470.110578
257	926520.0	9867.43693	1443687.6	2021-12-16	NA	NA	1453555.0	458.705254
258	926072.4	9632.54124	1444370.1	2021-12-17	NA	NA	1454002.6	447.571791
259	925635.7	9403.03124	1445036.3	2021-12-18	NA	NA	1454439.3	436.703950
260	925209.6	9178.79343	1445686.6	2021-12-19	NA	NA	1454865.4	426.095624
261	924793.9	8959.71620	1446321.4	2021-12-20	NA	NA	1455281.1	415.740833
262	924388.2	8745.68986	1446941.1	2021-12-21	NA	NA	1455686.8	405.633727
263	923992.5	8536.60659	1447545.9	2021-12-22	NA	NA	1456082.5	395.768576
264	923606.3	8332.36043	1448136.3	2021-12-23	NA	NA	1456468.7	386.139776
265	923229.6	8132.84731	1448712.6	2021-12-24	NA	NA	1456845.4	376.741842
266	922862.0	7937.96500	1449275.0	2021-12-25	NA	NA	1457213.0	367.569408
267	922503.4	7747.61313	1449824.0	2021-12-26	NA	NA	1457571.6	358.617223
268	922153.5	7561.69313	1450359.8	2021-12-27	NA	NA	1457921.5	349.880151
269	921812.2	7380.10826	1450882.7	2021-12-28	NA	NA	1458262.8	341.353168
270	921479.1	7202.76357	1451393.1	2021-12-29	NA	NA	1458595.9	333.031361
271	921154.2	7029.56588	1451891.2	2021-12-30	NA	NA	1458920.8	324.909922
272	920837.2	6860.42379	1452377.3	2021-12-31	NA	NA	1459237.8	316.984153
273	920528.0	6695.24762	1452851.8	2022-01-01	NA	NA	1459547.0	309.249456
274	920226.3	6533.94943	1453314.8	2022-01-02	NA	NA	1459848.7	301.701339
275	919931.9	6376.44298	1453766.6	2022-01-03	NA	NA	1460143.1	294.335407
276	919644.8	6222.64369	1454207.6	2022-01-04	NA	NA	1460430.2	287.147364
277	919364.7	6072.46867	1454637.9	2022-01-05	NA	NA	1460710.3	280.133011
278	919091.4	5925.83666	1455057.8	2022-01-06	NA	NA	1460983.6	273.288242
279	918824.8	5782.66801	1455467.6	2022-01-07	NA	NA	1461250.2	266.609044
280	918564.7	5642.88468	1455867.4	2022-01-08	NA	NA	1461510.3	260.091495
281	918310.9	5506.41021	1456257.6	2022-01-09	NA	NA	1461764.1	253.731762
282	918063.4	5373.16966	1456638.4	2022-01-10	NA	NA	1462011.6	247.526098
283	917822.0	5243.08965	1457010.0	2022-01-11	NA	NA	1462253.0	241.470841
284	917586.4	5116.09830	1457372.5	2022-01-12	NA	NA	1462488.6	235.562415
285	917356.6	4992.12520	1457726.3	2022-01-13	NA	NA	1462718.4	229.797322
286	917132.4	4871.10142	1458071.5	2022-01-14	NA	NA	1462942.6	224.172146
287	916913.7	4752.95946	1458408.3	2022-01-15	NA	NA	1463161.3	218.683551
288	916700.4	4637.63322	1458737.0	2022-01-16	NA	NA	1463374.6	213.328274
289	916492.3	4525.05802	1459057.6	2022-01-17	NA	NA	1463582.7	208.103130
290	916289.3	4415.17052	1459370.5	2022-01-18	NA	NA	1463785.7	203.005007
291	916091.3	4307.90875	1459675.8	2022-01-19	NA	NA	1463983.7	198.030863
292	915898.1	4203.21204	1459973.7	2022-01-20	NA	NA	1464176.9	193.177730
293	915709.6	4101.02104	1460264.3	2022-01-21	NA	NA	1464365.4	188.442706
294	915525.8	4001.27765	1460547.9	2022-01-22	NA	NA	1464549.2	183.822958
295	915346.5	3903.92506	1460824.6	2022-01-23	NA	NA	1464728.5	179.315717
296	915171.6	3808.90765	1461094.5	2022-01-24	NA	NA	1464903.4	174.918281
297	915001.0	3716.17106	1461357.9	2022-01-25	NA	NA	1465074.0	170.628009
298	914834.5	3625.66206	1461614.8	2022-01-26	NA	NA	1465240.5	166.442323
299	914672.2	3537.32863	1461865.5	2022-01-27	NA	NA	1465402.8	162.358706
300	914513.8	3451.11988	1462110.1	2022-01-28	NA	NA	1465561.2	158.374698
301	914359.3	3366.98603	1462348.7	2022-01-29	NA	NA	1465715.7	154.487898
302	914208.6	3284.87842	1462581.5	2022-01-30	NA	NA	1465866.4	150.695962
303	914061.6	3204.74945	1462808.6	2022-01-31	NA	NA	1466013.4	146.996600
304	913918.2	3126.55259	1463030.2	2022-02-01	NA	NA	1466156.8	143.387576
305	913778.4	3050.24235	1463246.4	2022-02-02	NA	NA	1466296.6	139.866708
306	913641.9	2975.77425	1463457.3	2022-02-03	NA	NA	1466433.1	136.431866
307	913508.8	2903.10482	1463663.1	2022-02-04	NA	NA	1466566.2	133.080967
308	913379.0	2832.19155	1463863.8	2022-02-05	NA	NA	1466696.0	129.811981
309	913252.4	2762.99290	1464059.6	2022-02-06	NA	NA	1466822.6	126.622924
310	913128.9	2695.46827	1464250.6	2022-02-07	NA	NA	1466946.1	123.511860
311	913008.4	2629.57797	1464437.0	2022-02-08	NA	NA	1467066.6	120.476900
312	912890.9	2565.28322	1464618.8	2022-02-09	NA	NA	1467184.1	117.516197
313	912776.3	2502.54612	1464796.2	2022-02-10	NA	NA	1467298.7	114.627950
314	912664.5	2441.32963	1464969.2	2022-02-11	NA	NA	1467410.5	111.810402
315	912555.4	2381.59757	1465138.0	2022-02-12	NA	NA	1467519.6	109.061835
316	912449.0	2323.31457	1465302.7	2022-02-13	NA	NA	1467626.0	106.380575
317	912345.3	2266.44609	1465463.3	2022-02-14	NA	NA	1467729.7	103.764985
318	912244.0	2210.95836	1465620.0	2022-02-15	NA	NA	1467831.0	101.213470
319	912145.3	2156.81841	1465772.9	2022-02-16	NA	NA	1467929.7	98.724472
320	912049.0	2103.99403	1465922.0	2022-02-17	NA	NA	1468026.0	96.296469
321	911955.1	2052.45373	1466067.5	2022-02-18	NA	NA	1468119.9	93.927977
322	911863.5	2002.16678	1466209.4	2022-02-19	NA	NA	1468211.5	91.617547
323	911774.1	1953.10315	1466347.8	2022-02-20	NA	NA	1468300.9	89.363767
324	911686.9	1905.23350	1466482.8	2022-02-21	NA	NA	1468388.1	87.165255
325	911601.9	1858.52919	1466614.5	2022-02-22	NA	NA	1468473.1	85.020664
326	911519.0	1812.96223	1466743.0	2022-02-23	NA	NA	1468556.0	82.928681
327	911438.1	1768.50528	1466868.4	2022-02-24	NA	NA	1468636.9	80.888023
328	911359.2	1725.13167	1466990.7	2022-02-25	NA	NA	1468715.8	78.897438
329	911282.3	1682.81532	1467109.9	2022-02-26	NA	NA	1468792.7	76.955704
330	911207.2	1641.53078	1467226.3	2022-02-27	NA	NA	1468867.8	75.061628
331	911134.0	1601.25319	1467339.8	2022-02-28	NA	NA	1468941.0	73.214049
332	911062.6	1561.95828	1467450.5	2022-03-01	NA	NA	1469012.4	71.411829
333	910992.9	1523.62235	1467558.5	2022-03-02	NA	NA	1469082.1	69.653861
334	910925.0	1486.22224	1467663.8	2022-03-03	NA	NA	1469150.0	67.939064
335	910858.7	1449.73537	1467766.6	2022-03-04	NA	NA	1469216.3	66.266384
336	910794.1	1414.13967	1467866.8	2022-03-05	NA	NA	1469280.9	64.634789
337	910731.0	1379.41359	1467964.6	2022-03-06	NA	NA	1469344.0	63.043277
338	910669.5	1345.53611	1468059.9	2022-03-07	NA	NA	1469405.5	61.490867
339	910609.6	1312.48669	1468153.0	2022-03-08	NA	NA	1469465.4	59.976601
340	910551.1	1280.24527	1468243.7	2022-03-09	NA	NA	1469523.9	58.499548
341	910494.0	1248.79230	1468332.2	2022-03-10	NA	NA	1469581.0	57.058796
342	910438.4	1218.10865	1468418.5	2022-03-11	NA	NA	1469636.6	55.653457
343	910384.1	1188.17569	1468502.8	2022-03-12	NA	NA	1469690.9	54.282664
344	910331.1	1158.97520	1468584.9	2022-03-13	NA	NA	1469743.9	52.945571
345	910279.5	1130.48941	1468665.0	2022-03-14	NA	NA	1469795.5	51.641352
346	910229.1	1102.70097	1468743.2	2022-03-15	NA	NA	1469845.9	50.369203
347	910180.0	1075.59294	1468819.4	2022-03-16	NA	NA	1469895.0	49.128337
348	910132.1	1049.14879	1468893.8	2022-03-17	NA	NA	1469942.9	47.917988
349	910085.3	1023.35239	1468966.3	2022-03-18	NA	NA	1469989.7	46.737408
350	910039.7	998.18798	1469037.1	2022-03-19	NA	NA	1470035.3	45.585867
351	909995.3	973.64020	1469106.1	2022-03-20	NA	NA	1470079.7	44.462654
352	909951.9	949.69404	1469173.4	2022-03-21	NA	NA	1470123.1	43.367073
353	909909.6	926.33486	1469239.1	2022-03-22	NA	NA	1470165.4	42.298447
354	909868.4	903.54836	1469303.1	2022-03-23	NA	NA	1470206.6	41.256115
355	909828.1	881.32060	1469365.6	2022-03-24	NA	NA	1470246.9	40.239431
356	909788.9	859.63796	1469426.5	2022-03-25	NA	NA	1470286.1	39.247767
357	909750.6	838.48715	1469485.9	2022-03-26	NA	NA	1470324.4	38.280508
358	909713.3	817.85521	1469543.9	2022-03-27	NA	NA	1470361.7	37.337055
359	909676.8	797.72949	1469600.4	2022-03-28	NA	NA	1470398.2	36.416825
360	909641.3	778.09763	1469655.6	2022-03-29	NA	NA	1470433.7	35.519246
361	909606.7	758.94759	1469709.4	2022-03-30	NA	NA	1470468.3	34.643763
362	909572.9	740.26760	1469761.8	2022-03-31	NA	NA	1470502.1	33.789833
363	909539.9	722.04618	1469813.0	2022-04-01	NA	NA	1470535.1	32.956927
364	909507.8	704.27215	1469862.9	2022-04-02	NA	NA	1470567.2	32.144528
365	909476.4	686.93456	1469911.6	2022-04-03	NA	NA	1470598.6	31.352132
366	909445.9	670.02276	1469959.1	2022-04-04	NA	NA	1470629.1	30.579249
367	909416.0	653.52633	1470005.4	2022-04-05	NA	NA	1470659.0	29.825399
368	909386.9	637.43513	1470050.6	2022-04-06	NA	NA	1470688.1	29.090113
369	909358.6	621.73924	1470094.7	2022-04-07	NA	NA	1470716.4	28.372936
370	909330.9	606.42899	1470137.7	2022-04-08	NA	NA	1470744.1	27.673423
371	909303.9	591.49496	1470179.6	2022-04-09	NA	NA	1470771.1	26.991139
372	909277.6	576.92794	1470220.5	2022-04-10	NA	NA	1470797.4	26.325661
373	909251.9	562.71893	1470260.4	2022-04-11	NA	NA	1470823.1	25.676576
374	909226.9	548.85919	1470299.3	2022-04-12	NA	NA	1470848.1	25.043480
375	909202.4	535.34016	1470337.2	2022-04-13	NA	NA	1470872.6	24.425981
376	909178.6	522.15349	1470374.2	2022-04-14	NA	NA	1470896.4	23.823694
377	909155.4	509.29105	1470410.3	2022-04-15	NA	NA	1470919.6	23.236247
378	909132.7	496.74488	1470445.6	2022-04-16	NA	NA	1470942.3	22.663273
379	909110.6	484.50725	1470479.9	2022-04-17	NA	NA	1470964.4	22.104416
380	909089.0	472.57059	1470513.4	2022-04-18	NA	NA	1470986.0	21.559330
381	909068.0	460.92753	1470546.1	2022-04-19	NA	NA	1471007.0	21.027676
382	909047.5	449.57085	1470577.9	2022-04-20	NA	NA	1471027.5	20.509123
383	909027.5	438.49356	1470609.0	2022-04-21	NA	NA	1471047.5	20.003348
384	909008.0	427.68878	1470639.3	2022-04-22	NA	NA	1471067.0	19.510038
385	908989.0	417.14985	1470668.9	2022-04-23	NA	NA	1471086.0	19.028885
386	908970.4	406.87023	1470697.7	2022-04-24	NA	NA	1471104.6	18.559590
387	908952.3	396.84357	1470725.9	2022-04-25	NA	NA	1471122.7	18.101862
388	908934.6	387.06365	1470753.3	2022-04-26	NA	NA	1471140.4	17.655416
389	908917.4	377.52443	1470780.1	2022-04-27	NA	NA	1471157.6	17.219973
390	908900.6	368.21999	1470806.2	2022-04-28	NA	NA	1471174.4	16.795264
391	908884.2	359.14457	1470831.6	2022-04-29	NA	NA	1471190.8	16.381024
392	908868.3	350.29256	1470856.4	2022-04-30	NA	NA	1471206.7	15.976994
393	908852.7	341.65845	1470880.7	2022-05-01	NA	NA	1471222.3	15.582925
394	908837.5	333.23691	1470904.3	2022-05-02	NA	NA	1471237.5	15.198569
395	908822.7	325.02270	1470927.3	2022-05-03	NA	NA	1471252.3	14.823689
396	908808.2	317.01075	1470949.8	2022-05-04	NA	NA	1471266.8	14.458051
397	908794.1	309.19607	1470971.7	2022-05-05	NA	NA	1471280.9	14.101427
398	908780.4	301.57383	1470993.1	2022-05-06	NA	NA	1471294.6	13.753595
399	908766.9	294.13929	1471013.9	2022-05-07	NA	NA	1471308.1	13.414339
400	908753.9	286.88784	1471034.3	2022-05-08	NA	NA	1471321.1	13.083448
401	908741.1	279.81499	1471054.1	2022-05-09	NA	NA	1471333.9	12.760714
402	908728.6	272.91633	1471073.4	2022-05-10	NA	NA	1471346.4	12.445939
403	908716.5	266.18760	1471092.3	2022-05-11	NA	NA	1471358.5	12.138924
404	908704.7	259.62461	1471110.7	2022-05-12	NA	NA	1471370.3	11.839480
405	908693.1	253.22328	1471128.7	2022-05-13	NA	NA	1471381.9	11.547420
406	908681.9	246.97965	1471146.2	2022-05-14	NA	NA	1471393.1	11.262561
407	908670.9	240.88983	1471163.2	2022-05-15	NA	NA	1471404.1	10.984727
408	908660.2	234.95004	1471179.9	2022-05-16	NA	NA	1471414.8	10.713744
409	908649.7	229.15660	1471196.1	2022-05-17	NA	NA	1471425.3	10.449443
410	908639.5	223.50589	1471212.0	2022-05-18	NA	NA	1471435.5	10.191660
411	908629.6	217.99442	1471227.4	2022-05-19	NA	NA	1471445.4	9.940235
412	908619.9	212.61875	1471242.5	2022-05-20	NA	NA	1471455.1	9.695009
413	908610.4	207.37554	1471257.2	2022-05-21	NA	NA	1471464.6	9.455832
414	908601.2	202.26154	1471271.5	2022-05-22	NA	NA	1471473.8	9.222553
415	908592.2	197.27357	1471285.5	2022-05-23	NA	NA	1471482.8	8.995027
416	908583.4	192.40852	1471299.2	2022-05-24	NA	NA	1471491.6	8.773113
417	908574.9	187.66336	1471312.5	2022-05-25	NA	NA	1471500.1	8.556672
418	908566.5	183.03516	1471325.4	2022-05-26	NA	NA	1471508.5	8.345569
419	908558.4	178.52102	1471338.1	2022-05-27	NA	NA	1471516.6	8.139673
420	908550.5	174.11815	1471350.4	2022-05-28	NA	NA	1471524.5	7.938855
421	908542.7	169.82379	1471362.5	2022-05-29	NA	NA	1471532.3	7.742990
422	908535.2	165.63529	1471374.2	2022-05-30	NA	NA	1471539.8	7.551956
423	908527.8	161.55003	1471385.7	2022-05-31	NA	NA	1471547.2	7.365634
424	908520.6	157.56548	1471396.8	2022-06-01	NA	NA	1471554.4	7.183908
425	908513.6	153.67915	1471407.7	2022-06-02	NA	NA	1471561.4	7.006665
426	908506.8	149.88862	1471418.3	2022-06-03	NA	NA	1471568.2	6.833793
427	908500.1	146.19154	1471428.7	2022-06-04	NA	NA	1471574.9	6.665185
428	908493.6	142.58560	1471438.8	2022-06-05	NA	NA	1471581.4	6.500737
429	908487.3	139.06856	1471448.7	2022-06-06	NA	NA	1471587.7	6.340345
430	908481.1	135.63823	1471458.3	2022-06-07	NA	NA	1471593.9	6.183909
431	908475.1	132.29248	1471467.7	2022-06-08	NA	NA	1471599.9	6.031333
432	908469.2	129.02921	1471476.8	2022-06-09	NA	NA	1471605.8	5.882520
433	908463.4	125.84641	1471485.7	2022-06-10	NA	NA	1471611.6	5.737378
434	908457.8	122.74208	1471494.4	2022-06-11	NA	NA	1471617.2	5.595817
435	908452.4	119.71430	1471502.9	2022-06-12	NA	NA	1471622.6	5.457747
436	908447.1	116.76117	1471511.2	2022-06-13	NA	NA	1471627.9	5.323084
437	908441.9	113.88086	1471519.3	2022-06-14	NA	NA	1471633.1	5.191743
438	908436.8	111.07158	1471527.1	2022-06-15	NA	NA	1471638.2	5.063642
439	908431.9	108.33157	1471534.8	2022-06-16	NA	NA	1471643.1	4.938701
440	908427.0	105.65912	1471542.3	2022-06-17	NA	NA	1471648.0	4.816842
441	908422.3	103.05258	1471549.6	2022-06-18	NA	NA	1471652.7	4.697990
442	908417.8	100.51032	1471556.7	2022-06-19	NA	NA	1471657.2	4.582070
443	908413.3	98.03075	1471563.7	2022-06-20	NA	NA	1471661.7	4.469009
444	908408.9	95.61234	1471570.5	2022-06-21	NA	NA	1471666.1	4.358738
445	908404.7	93.25356	1471577.1	2022-06-22	NA	NA	1471670.3	4.251188
446	908400.5	90.95296	1471583.5	2022-06-23	NA	NA	1471674.5	4.146290
447	908396.5	88.70909	1471589.8	2022-06-24	NA	NA	1471678.5	4.043981
448	908392.5	86.52057	1471595.9	2022-06-25	NA	NA	1471682.5	3.944196
449	908388.7	84.38602	1471601.9	2022-06-26	NA	NA	1471686.3	3.846873
450	908385.0	82.30412	1471607.7	2022-06-27	NA	NA	1471690.0	3.751950
451	908381.3	80.27357	1471613.4	2022-06-28	NA	NA	1471693.7	3.659370
452	908377.7	78.29310	1471619.0	2022-06-29	NA	NA	1471697.3	3.569074
453	908374.2	76.36148	1471624.4	2022-06-30	NA	NA	1471700.8	3.481006