Snippets
Scenarios
Random draws
Simulate random draws
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# This creates a number of
# scenarios by drawing random
# values from four different
# distributions. Afterwards the
# the results are saved to a
# scenarios document item
# named "Random"
# Lower precision a bit to
# speed up simulations
accuracy( 20 );
# Number of scenarios
numS = 10;
# Number of draws per scenario
numD = 50;
# Date offset
# This can be 0 for indexed
# (no frequency) values.
# Otherwise use the first date of
# the simulation sample here.
I = 0;
# For non-index values, uncomment
# the line below
# I = dateIndex( 2017, 1, 1 ) - 1;
# Observation matrices
normObs = matrix( numD,
1 + numS );
chiSObs = matrix( numD,
1 + numS );
fObs = matrix( numD,
1 + numS );
betaObs = matrix( numD,
1 + numS );
s = 1;
while ( s <= numS,
i = 1,
while( i <= numD,
# Set date column value
# Note that we use the
# date offset to correct
if ( s == 1,
{ normObs[ i;1 ] = i + I,
chiSObs[ i;1 ] = i + I,
fObs [ i;1 ] = i + I,
betaObs[ i;1 ] = i + I } ),
# Simulate
normObs[ i;1+s ] =
shockAdd( 0, "norm", 0, 1 ),
chiSObs[ i;1+s ] =
shockAdd( 0, "chisq", 1, 0 ),
fObs [ i;1+s ] =
shockAdd( 0, "f", 5, 7 ),
betaObs[ i;1+s ] =
shockAdd( 0, "beta", 2, 2 ),
i = i + 1 ),
# Progress
calcStatus( "Scenario " + s + "/" + numS ),
s = s + 1 );
# Now create new scenarios
scenariosAdd( "Random",
"Normal", normObs,
"Chi square", chiSObs,
"F", fObs,
"Beta", betaObs );
Random walks
Simulate random walks
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# This creates a number of
# scenarios by simulating
# random walks using four
# different distributions.
# The results are saved to a
# scenarios document item
# named "Random walks"
# Lower precision a bit to
# speed up simulations
accuracy( 20 );
# Number of scenarios
numS = 20;
# Number of draws per scenario
numD = 100;
# Date offset
# This can be 0 for indexed
# (no frequency) values.
# Otherwise use the first date of
# the simulation sample here.
I = 0;
# For non-index values, uncomment
# the line below
# I = dateIndex( 2017, 1, 1 ) - 1;
# Observation matrices
normObs = matrix( numD,
1 + numS );
chiSObs = matrix( numD,
1 + numS );
fObs = matrix( numD,
1 + numS );
betaObs = matrix( numD,
1 + numS );
s = 1;
while ( s <= numS,
# First observation is 0
# and random walk starting
# point for all walks
i = 2,
while( i <= numD,
# Set date column value
# Note that we use the
# date offset to correct
if ( s == 1,
{ normObs[ i;1 ] = i + I,
chiSObs[ i;1 ] = i + I,
fObs [ i;1 ] = i + I,
betaObs[ i;1 ] = i + I } ),
# Simulate
normObs[ i;1+s ] =
shockAdd( normObs[ i-1;1+s ],
"norm", 1, 1 ),
chiSObs[ i;1+s ] =
shockAdd( chiSObs[ i-1;1+s ],
"chisq", 1, 0 ),
fObs [ i;1+s ] =
shockAdd( fObs [ i-1;1+s ],
"f", 5, 7 ),
betaObs[ i;1+s ] =
shockAdd( betaObs[ i-1;1+s ],
"beta", 2, 2 ),
i = i + 1 ),
# Progress
calcStatus( "Scenario " + s + "/" + numS ),
s = s + 1 );
# We still need to set the date
# for the first observation. Since
# all random walks started at 0
# we started the loop above at 2
# and not 1, so as a side effect
# the first date has not been
# loaded yet. This only matters
# when I, the date offset, contains
# a non-zero value.
normObs[ 1;1 ] = I + 1;
chiSObs[ 1;1 ] = I + 1;
fObs [ 1;1 ] = I + 1;
betaObs[ 1;1 ] = I + 1;
scenariosAdd( "Random walks",
"Normal", normObs,
"Chi square", chiSObs,
"F", fObs,
"Beta", betaObs );
Virus infections
Simulate pandemic (pareto)
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# Virus infections simulation
# Use pareto to simulate contact
# rate per week
# Lower precision a bit to
# speed up simulations
accuracy( 20 );
# Numbers get big so change format
format( 3, 30 );
# We perform two simulations
# First we do the baseline, then we
# simulate scenarios
# So we loop the same code twice
loop = 0;
while( loop < 2,
{
# Number of scenarios
# If this is 1, the baseline will be simulated
# Otherwise this holds the number of scenarios
# to simulate
# Set to 1 if loop is 0 and 100 if loop is not
if( loop, numS = 100, numS = 1 );
# Weeks per scenario
numT = 52;
# Starting date
t0 = dateIndex ( 2020, 1, 1 );
# Number of days per interval
td = 7;
# Assumptions - fixed
p0 = 100000000; # Initial population 100m
i0 = 10; # Initial infections
pRec = 0.9925; # Probability of recovery (cat III)
tIll = 3; # Weeks to recover
cRBase = 1; # Baseline contact rate per week
# Assumptions - variable
# Pareto scale and shape to simulate
# contact rate
# The shape 1.16 is standard 80-20 for pareto
pAvg = 1; # the average is 1
pShape = 1.16; # Arbitrary, mean is cNum
pScale = pAvg * ( pShape - 1 ) / pShape;
pFast = 1 / pShape; # Speed up calculations later
# Variables below
# tPop - total population
# tInf - total infected
# nInf - newly infected
# tSuc - total susceptible
# nRec - newly recovered
# tRec - total recovered
# nDec - newly deceased
# tDec - total deceased
# Matrices
# First column is date
nInfMat = matrix( numT, 1 + numS );
tRecMat = matrix( numT, 1 + numS );
tDecMat = matrix( numT, 1 + numS );
cRateMat = matrix( numT, 1 + numS );
# Simulate
s = 1; # Scenario
while ( s <= numS,
{
# Initialise for this scenario
tPop = p0; # Population
tInf = i0; # Total infected
tRec = 0; # Total recovered
tDec = 0; # Total deceased
nRec = 0; # Newly recovered
nDec = 0; # Newly deceased
T = t0; # Date index
t = 1; # Counter
while ( t <= numT,
{
# Population
tPop = tPop - nDec;
if( numS > 1,
# This is not baseline
# Simulate contact rate
cR = pScale / ( 1 - randFast () ) ^ pShape,
# Else this is baseline
# Use constant contact rate
cR = cRBase
);
# Susceptible population
tSuc = tPop - tInf - tRec;
# Newly infected
nInf = 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
nit = nInfMat[ t - tIll; s + 1 ];
# Recovered
nRec = nit * pRec;
# Deceased
nDec = nit - nRec;
} );
# Total infected
tInf = tInf + nInf - nRec - nDec;
# Total recovered and deceassed
tRec = tRec + nRec;
tDec = tDec + nDec;
# Load date column
if ( s == 1,
{
nInfMat [ 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;
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
if( numS > 1,
calcStatus( "Scenario " + s + "/" + numS ),
calcStatus( "Baseline" ) );
} );
# Simulation name
if( numS > 1,
# then
name = "Pandemic scenarios",
# else
name = "Pandemic baseline" );
# Prevent clutter
scenariosDelete( name );
# Create scenarios
scenariosAdd(
name,
"New Infections", nInfMat,
"Recovered", tRecMat,
"Deceased", tDecMat,
"Contact Rate", cRateMat );
# Print some baseline variables
if ( numS == 1, (
print( "Baseline - new infections" );
print( matrixAsSeries( nInfMat ) );
print( "Baseline - recovered" );
print( matrixAsSeries( tRecMat ) );
print( "Baseline - deceased" );
print( matrixAsSeries( tDecMat ) );
print( "Baseline - contact rate" );
print( matrixAsSeries( cRateMat ) );
# Also save for later use
dataReplace( "NewInf0", matrixAsSeries( nInfMat ) );
dataReplace( "Recovered0", matrixAsSeries( tRecMat ) );
dataReplace( "Deceased0", matrixAsSeries( tDecMat ) );
dataReplace( "Contact0", matrixAsSeries( cRateMat ) );
) );
# Next loop
loop = loop + 1
} );
Fast Pareto only
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# Use the Pareto distribution to
# simulated the contact rate.
# Drawing random numbers from
# pareto can be sped up.
# Lower precision a bit to
# speed up simulations
accuracy( 20 );
# Number of scenarios
numS = 1000;
# Number of draws per scenario
numD = 5;
# Date offset
# We use index based observations
# so any integer will do
I = 0;
# Observation matrices
# Fist column is for date / index
paretoObs = matrix( numD, 1 + numS );
# Some parameters we need for
# the fast method
# Pareto
pAvg = 1; # Average
pShape = 1.16;
pScale = pAvg * ( pShape - 1 ) / pShape;
pFast = 1 / pShape; # For later
# Start simulating
s = 1;
while ( s <= numS,
i = 1,
while( i <= numD,
# Set date column value
# Note that we use the
# date offset to correct
if ( s == 1,
paretoObs[ i;1 ] = i + I ),
# Simulate
paretoObs[ i;1+s ] = pScale / ( 1 - randFast () ) ^ pFast,
i = i + 1 ),
# Progress
calcStatus( "Scenario " + s + "/" + numS ),
s = s + 1 );
# Delete any existing simulation
# to prevent clutter
scenariosDelete( "ContactSimPareto" );
# Now create new scenarios
scenariosAdd( "ContactSimPareto",
"Pareto", paretoObs );
Virus infections
Simulate pandemic (p+b)
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# Virus infections simulation
# Using both beta and pareto
# distributions
# beta to simulate probability
# of transmitting disease
# pareto to simulate number of
# contacts per week
# Together this simulates the
# contact rate.
# Lower precision a bit to
# speed up simulations
accuracy( 20 );
# Numbers get big so change format
format( 3, 30 );
# We perform two simulations
# First we do the baseline, then we
# simulate scenarios
# So we loop the same code twice
loop = 0;
while( loop < 2,
{
# Number of scenarios
# If this is 1, the baseline will be simulated
# Otherwise this holds the number of scenarios
# to simulate
# Set to 1 if loop is 0 and 100 if loop is not
if( loop, numS = 100, numS = 1 );
# Weeks per scenario
numT = 52;
# Starting date
t0 = dateIndex ( 2020, 1, 1 );
# Number of days per interval
td = 7;
# Assumptions - fixed
p0 = 100000000; # Initial population 100m
i0 = 10; # Initial infections
pRec = 0.9925; # Probability of recovery (cat III)
tIll = 3; # Weeks to recover
# Assumptions - variable with these as mean
pInf = 0.2; # Probability of infecting a contact per interval
cNum = 5; # Number of people contacted per interval
cRBase = pInf * cNum; # Baseline contact rate
# Beta parameters used to simulate infection rate
pAlpha = 10 * pInf; # Arbitrary
pBeta = 10 - pAlpha; # Mean of distro is pInf
# If pAlpha is an integer, we can speed up the random number generation
bFast = round( pAlpha ) == pAlpha;
bCount = pAlpha + pBeta - 1;
bIndex = pAlpha;
# Pareto scale and shape to simulate contacts
# The shape 1.16 is standard 80-20 for pareto
# It is also a bit wild, make the shape bigger
# for a better behaving result
pShape = 1.16; # Arbitrary, mean is cNum
pScale = cNum * ( pShape - 1 ) / pShape;
pFast = 1 / pShape; # Speed up calculations later
# Variables below
# tPop - total population
# tInf - total infected
# nInf - newly infected
# tSuc - total susceptible
# nRec - newly recovered
# tRec - total recovered
# nDec - newly deceased
# tDec - total deceased
# Matrices
# First column is date
nInfMat = matrix( numT, 1 + numS );
tRecMat = matrix( numT, 1 + numS );
tDecMat = matrix( numT, 1 + numS );
cRateMat = matrix( numT, 1 + numS );
# Simulate
s = 1; # Scenario
while ( s <= numS,
{
# Initialise for this scenario
tPop = p0; # Population
tInf = i0; # Total infected
tRec = 0; # Total recovered
tDec = 0; # Total deceased
nRec = 0; # Newly recovered
nDec = 0; # Newly deceased
T = t0; # Date index
t = 1; # Counter
while ( t <= numT,
{
# Population
tPop = tPop - nDec;
if( numS > 1,
# This is not baseline
# Simulate random contacts
(
# Simulate contacts
nC = pScale / ( 1 - randFast () ) ^ pShape;
# Simulate infection rate
if( bFast,
# then
if ( bIndex == 1,
# even faster if a == 1
pI = min( randFast( 1, bCount ) ),
# fast way, a'th variable from sorted a+b-1 sample
(
uRnd = sort( randFast( 1, bCount ) ),
pI = uRnd[ bIndex ]
)
),
# else
# slow way, inverse beta
pI = betaQnt( randFast (), pAlpha, pBeta )
);
# Contact rate
cR = nC * pI;
),
# else this is baseline
(
cR = cRBase;
) );
# Susceptible population
tSuc = tPop - tInf - tRec;
# Newly infected
nInf = 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
nit = nInfMat[ t - tIll; s + 1 ];
# Recovered
nRec = nit * pRec;
# Deceased
nDec = nit - nRec;
} );
# Total infected
tInf = tInf + nInf - nRec - nDec;
# Total recovered and deceassed
tRec = tRec + nRec;
tDec = tDec + nDec;
# Load date column
if ( s == 1,
{
nInfMat [ 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;
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
if( numS > 1,
calcStatus( "Scenario " + s + "/" + numS ),
calcStatus( "Baseline" ) );
} );
# Simulation name
if( numS > 1,
# then
name = "Pandemic scenarios",
# else
name = "Pandemic baseline" );
# Prevent clutter
scenariosDelete( name );
# Create scenarios
scenariosAdd(
name,
"New Infections", nInfMat,
"Recovered", tRecMat,
"Deceased", tDecMat,
"Contact Rate", cRateMat );
# Print some baseline variables
if ( numS == 1, (
print( "Baseline - new infections" );
print( matrixAsSeries( nInfMat ) );
print( "Baseline - recovered" );
print( matrixAsSeries( tRecMat ) );
print( "Baseline - deceased" );
print( matrixAsSeries( tDecMat ) );
print( "Baseline - contact rate" );
print( matrixAsSeries( cRateMat ) );
# Also save for later use
dataReplace( "NewInf0", matrixAsSeries( nInfMat ) );
dataReplace( "Recovered0", matrixAsSeries( tRecMat ) );
dataReplace( "Deceased0", matrixAsSeries( tDecMat ) );
dataReplace( "Contact0", matrixAsSeries( cRateMat ) );
) );
# Next loop
loop = loop + 1
} );
Fast Pareto and Beta
# Copyright © Ferra Solutions (Pty) Ltd. All rights reserved.
# Using both beta and pareto
# distributions
# beta to simulate probability
# of transmitting disease
# pareto to simulate number of
# contacts per week
# Together this simulates the
# contact rate.
# Drawing random numbers from
# pareto and beta distributions
# can be sped up. This creates
# a distribution based on the
# fast method so that it can be
# verified and also to investigate
# the outcome of the simulation.
# Lower precision a bit to
# speed up simulations
accuracy( 20 );
# Number of scenarios
numS = 1000;
# Number of draws per scenario
numD = 5;
# Date offset
# We use index based observations
# so any integer will do
I = 0;
# Observation matrices
# Fist column is for date / index
paretoObs = matrix( numD,
1 + numS );
betaObs = matrix( numD,
1 + numS );
# Some parameters we need for
# the fast method
# Pareto
pAvg = 5; # Average
pShape = 1.16;
pScale = pAvg * ( pShape - 1 ) / pShape;
pFast = 1 / pShape; # For later
# Beta
bNum = 9;
bIdx = 2;
# Start simulating
s = 1;
while ( s <= numS,
i = 1,
while( i <= numD,
# Set date column value
# Note that we use the
# date offset to correct
if ( s == 1,
{ paretoObs[ i;1 ] = i + I,
betaObs [ i;1 ] = i + I } ),
# Simulate
paretoObs[ i;1+s ] = pScale / ( 1 - randFast () ) ^ pFast,
betaDraw = sort( randFast( 1, bNum ) ),
betaObs [ i;1+s ] = betaDraw[ bIdx ],
i = i + 1 ),
# Progress
calcStatus( "Scenario " + s + "/" + numS ),
s = s + 1 );
# Calculate contact rate
# Element by element multiply
cRateObs = matrixElMul( paretoObs, betaObs );
# Now fix up the dates
# Replace the first column with that of betaObs
cRateObs[ 1:numD;1 ] = betaObs[ 1:numD;1 ];
# Delete any existing simulation
# to prevent clutter
scenariosDelete( "ContactSim" );
# Now create new scenarios
scenariosAdd( "ContactSim",
"Pareto", paretoObs,
"Beta", betaObs,
"Contact Rate", cRateObs );