18a sir vs reality

In this video I discuss a number of sir refinements and calibrations to get the model to approximate reality. I also draw some conclusions - see the video and let me know in the comments if you agree!

The fdm simulation is given below. It improves on earlier versions by also discretising the population - the population is divided up into 100m actual units. My plan here was to simulate infections using a poisson binomial but as I also tried to incorporate the spatial interaction I went hardcore and rather programmed it from the ground up. So what you see in the video is somewhat disconnected from fdm.

# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.

# Lower precision a bit to
# speed up simulations
accuracy( 20 );

# Given a range [ l, h ] and a divisor m
# find the number of integers l, l+1, ..., h
# that are divisible by m
define( "countMultiples", "l", "h", "m",
  {
    m0 = ceil  ( l / m );
    m1 = floor ( h / m );
    m1 - m0 + 1;
  } );

# Weeks per scenario
numT   =   52;

# Starting date
t0     = dateIndex ( 2020, 1, 1 );

# Number of days per interval
td     = 7;

# Assumptions
p0     = 100000000;    # Initial population
i0     = 10;           # Initial infections
inf1   = 50;           # For every inf1 infected
                       # we will have 1 fatality
tIll   = 3;            # Weeks to recover
cRBase = 1;            # Baseline contact rate
numAlt = 100;          # Number of alternative
                       # scenarios

# Parato parameters
pShape = 1.16;         # Arbitrary 80:20
pFast  = 1 / pShape;   # Speed up calculations
                       # Scale based on average
pScale = cRBase * ( pShape - 1 ) / pShape;

# We perform two simulations

# First we do the baseline, then we
# simulate variable contact rate
# scenarios

# So we loop the same code twice
loop = 0;

while( loop < 2, {

  # Variables below

  # tPop - total population
  # tInf - total infected
  # cInf - cumulative infected
  # nInf - newly infected
  # tSuc - total susceptible
  # nRec - newly recovered
  # tRec - total recovered
  # nDec - newly deceased
  # tDec - total deceased

  # Number of scenarios
  numS = if ( loop, numAlt, 1 );

  # Matrices
  # First column is date
  nInfMat  = matrix( numT, 1 + numS );
  tInfMat  = matrix( numT, 1 + numS );
  cInfMat  = matrix( numT, 1 + numS );
  tRecMat  = matrix( numT, 1 + numS );
  tDecMat  = matrix( numT, 1 + numS );
  cRateMat = matrix( numT, 1 + numS );

  # Simulate
  s = 1; # Scenario number

  while ( s <= numS,
  {
    # Initialise for this scenario
    tPop = p0; # Population
    tInf =  0; # Total infected
    cInf =  0; # Cumulative infected
    tRec =  0; # Total recovered
    tDec =  0; # Total deceased
    nRec =  0; # Newly recovered
    nDec =  0; # Newly deceased

    T    = t0; # Date index
    t    =  1; # Time bucket

    while ( t <= numT,
    {
      # Population adjusted
      # for deceased
      tPop = tPop - nDec;

      if( numS > 1,

        # This is not baseline
        # Simulate random contacts
        cR = pScale / ( 1 - randFast () ) ^ pFast,

        # else this is baseline
        cR = cRBase
      );

      # Susceptible population
      tSuc = tPop - tInf - tRec;

      # Newly infected
      if ( t == 1,
        # First time it is our assumption
        nInf = i0,
	# Thereafter it is the rounded number
        nInf = round( cR * tSuc * tInf / tPop ) );

      # Can not exceed susceptible
      if ( nInf > tSuc, nInf = tSuc );

      # Recovered and deceased
      if ( t > tIll,
      # then
      {
        # Newly infected tIll time ago
        nia  = nInfMat[ t - tIll; s + 1 ];

        # Cumulative infected tIll time ago
	cia  = cInfMat[ t - tIll; s + 1 ];

        # Fatality range
	f0   = cia - nia + 1;
	f1   = f0  + nia - 1;

        # Fatalities / newly deceased
	nDec = call( "countMultiples", f0, f1, inf1 );

        # Recovered - infected less fatalities
        nRec = nia - nDec;

      } );

      # Total infected
      tInf = tInf + nInf - nRec - nDec;

      # Cumulative infected
      cInf = cInf + nInf;

      # Total recovered and deceased
      tRec = tRec + nRec;
      tDec = tDec + nDec;

      # Load date column
      if ( s == 1,
      {
        nInfMat [ t; s ] = T;
	tInfMat [ t; s ] = T;
	cInfMat [ t; s ] = T;
        tRecMat [ t; s ] = T;
        tDecMat [ t; s ] = T;
        cRateMat[ t; s ] = T;
      } );

      # Update matrices
      # First column is date
      nInfMat [ t; s + 1 ] = nInf;
      tInfMat [ t; s + 1 ] = tInf;
      cInfMat [ t; s + 1 ] = cInf;
      tRecMat [ t; s + 1 ] = tRec;
      tDecMat [ t; s + 1 ] = tDec;
      cRateMat[ t; s + 1 ] = cR;

      # Move to next time
      T = T + td;
      t = t + 1;
    } ),

    # Move to next scenario
    s = s + 1;

    # Progress
    format( 0 );
    calcStatus( "Scenario " + s + "/" + numS );
  } );

  # Simulation name
  if( numS > 1,
      # then
      name = "Scenarios",
      # else
      name = "Baseline" );

  # Prevent clutter
  scenariosDelete( name );

  # Create scenarios
  scenariosAdd(
    name,
      "New Infections", nInfMat,
      "Total Inf",      tInfMat,
      "Cumulative Inf", cInfMat,
      "Recovered",      tRecMat,
      "Deceased",       tDecMat,
      "Contact Rate",   cRateMat );

  # Print results
  if( loop == 0, {
      format( 2, 30 );
      print( "New infections" );
      print( matrixAsSeries( nInfMat ) );
      print( "Total infections" );
      print( matrixAsSeries( tInfMat ) );
      print( "Cumulative infections" );
      print( matrixAsSeries( cInfMat ) );
      print( "Total recovered" );
      print( matrixAsSeries( tRecMat ) );
      print( "Total deceased" );
      print( matrixAsSeries( tDecMat ) );
  } );

  # Next loop
  loop = loop + 1 } );