Simulation and Modeling

Stochastic activities

A stochastic activity is an activity whose duration or outcome is not fixed but random, described by a probability distribution. The time to complete it cannot be predicted exactly, e.g. the service time of a customer at a bank or the arrival of a call. Such activities are simulated by sampling values from distributions (exponential, normal, uniform, Poisson) using random numbers. In contrast, a deterministic activity has a fixed, known duration.

Principles of modelling

1. Block building: the model should be built from well-defined, independent blocks (subsystems), each with clearly defined inputs and outputs, so that blocks can be tested and replaced separately.
2. Relevance: include only the information relevant to the questions being studied; irrelevant detail adds complexity without value.
3. Accuracy: the accuracy of information included must be consistent with the purpose and with the quality of available data.
4. Aggregation: group parts of the real system into aggregate units where the detail is unnecessary, to keep the model manageable.
5. Fallibility: a model is only an approximation; it cannot be perfect and must be validated against the real system.
6. Simplicity (parsimony): start simple and add detail gradually; prefer the simplest model that meets the objective.
7. Iteration: modelling is an evolving process of build, test, validate and refine.
8. Verification and validation: check that the model is implemented correctly and that it represents the real system.
9. Clear objectives: define the purpose, boundary, variables, parameters and performance measures before building.

Distributed lag model

A distributed lag model is a dynamic model in which the effect of an explanatory variable on the dependent variable is spread over several time periods, so the current value of $Y$ depends on the present and past values of $X$:

$Y_t = \alpha + \beta_0 X_t + \beta_1 X_{t-1} + \beta_2 X_{t-2} + \dots + \beta_k X_{t-k} + u_t$

Examples: the effect of advertising or of income on consumption, or price on supply, takes several periods. The coefficients $\beta_i$ show the effect after $i$ periods, and their sum is the long-run multiplier. Special forms are the Koyck (geometric) lag and the Almon polynomial lag.

Monte Carlo estimate of an integral

The integral printed in the paper is not reproduced here, so the method is shown on a typical example, $I=\int_0^2 x^2\,dx$ (exact value $8/3 = 2.667$).

Method: for $I=\int_a^b f(x)\,dx$, generate uniform random numbers $R_i$ in (0,1), set $x_i = a + (b-a)R_i$ and estimate

$I \approx (b-a)\,\frac{1}{N}\sum_{i=1}^{N} f(x_i)$

Take $N=5$ random numbers 0.15, 0.62, 0.37, 0.89, 0.51 with $a=0$, $b=2$:

Sum = 6.384, mean = 1.2768, so

$I \approx 2 \times 1.2768 = 2.5536$

This is close to the exact 2.667, and the estimate improves as $N$ increases (the error falls as $1/\sqrt{N}$). The same steps are applied to the given integral by changing $a$, $b$ and $f(x)$.

Distributed lag model

A distributed lag model expresses a dependent variable in terms of the current and past values of explanatory variables, because the response to a change is spread over time:

$Y_t = \alpha + \beta_0X_t + \beta_1X_{t-1} + \dots + \beta_kX_{t-k} + u_t$

Market model with lag (cobweb)

Producers decide supply on the price of the previous period, since production takes time, while demand responds to the current price:

$\text{Demand: } Q_t^d = a - bP_t \qquad \text{Supply: } Q_t^s = c + dP_{t-1}$

In equilibrium $Q^d_t=Q^s_t$:

$a - bP_t = c + dP_{t-1} \;\Rightarrow\; P_t = \frac{a-c}{b} - \frac{d}{b}P_{t-1}$

With $a=12$, $b=2$, $c=5$, $d=0.9$:

$P_t = 3.5 - 0.45\,P_{t-1}$

Equilibrium price: $P^* = \dfrac{a-c}{b+d} = \dfrac{7}{2.9} = 2.414$.

Starting from $P_0 = 5$:

Cobweb graph: plot price on the vertical axis and quantity on the horizontal. Draw the demand line $P = 6 - 0.5Q$ and the supply line $P=(Q-5)/0.9$. From $P_0=5$ move horizontally to the supply curve (quantity supplied next period), then vertically to the demand curve to read the new price $P_1=1.25$, then back to supply, and so on. The path spirals inward towards the intersection at $P^*\approx2.41$.

Conclusion: since the slope ratio $d/b = 0.45 < 1$ (the coefficient of $P_{t-1}$ has magnitude below 1), the price oscillates with decreasing amplitude and converges: the market is stable (convergent cobweb).

Analog computer

An analog computer represents variables of a problem by continuously varying physical quantities (usually voltages) and solves differential equations by electrical circuits that behave like the system modelled. It works in real time, in parallel and continuously, as opposed to a digital computer, which works on discrete numbers.

Components

1. Summer (adder): produces the sum of input voltages (with sign inversion), $e_o = -(e_1+e_2+\dots)$.
2. Integrator: gives the time integral of its input, $e_o = -\frac{1}{RC}\int e_i\,dt$; the key element for differential equations.
3. Coefficient potentiometer: multiplies a signal by a constant between 0 and 1.
4. Inverter (sign changer): reverses the sign.
5. Multiplier and function generators: form products and nonlinear functions.
6. Operational amplifiers: high-gain amplifiers used to build the above.
7. Comparators, switches and initial-condition circuits.
8. Output devices: oscilloscope, plotter or recorder.

Block diagram for the liver

The liver is modelled with a compartment model, e.g. a drug (or substance) in the blood and liver. Let $x$ = amount in blood and $y$ = amount in liver, with rate constants $k_1$ (blood to liver) and $k_2$ (liver elimination), and an input $u(t)$:

$\frac{dx}{dt} = u - k_1 x, \qquad \frac{dy}{dt} = k_1x - k_2y$

Analog setup:

1. Input $u$ and the feedback $-k_1x$ go into a summer; its output is $dx/dt$.
2. An integrator (initial value $x_0$) produces $-x$; a potentiometer set to $k_1$ forms $k_1x$, which is fed back to the first summer through an inverter, and also forwarded to the second stage.
3. Second stage: $k_1x$ and $-k_2y$ feed a summer giving $dy/dt$, followed by an integrator (initial value $y_0$) giving $y$; a potentiometer set to $k_2$ returns $k_2y$ as feedback.
4. The outputs $x(t)$ and $y(t)$ are recorded on a plotter.

Block flow: $u \to$ [Summer] $\to$ [Integrator] $\to x \to$ [Pot $k_1$] $\to$ feedback to the summer and forward to [Summer] $\to$ [Integrator] $\to y \to$ [Pot $k_2$] $\to$ feedback.

CSMP

CSMP (Continuous System Modelling Program) is a simulation language developed by IBM for modelling continuous systems described by differential equations. It was designed to replace analog computers: the user writes the equations in block-diagram-like statements, the processor sorts them into the correct computational order and integrates them numerically, and prints or plots the results. Its statements are order-independent (the program sequences them).

Statement kinds

1. Structure statements: describe the model using functional blocks such as INTGRL, and algebraic expressions.
   • Y = INTGRL(Y0, DY) integrates DY with initial value Y0.
   • DY = -K*Y + U is an equation for a derivative.
   • Function blocks: STEP(T0), RAMP(T0), IMPULS, DELAY, LIMIT, FCNSW.
2. Data statements: assign values to parameters, constants and initial conditions.
   • PARAMETER K = 0.5, CONSTANT Y0 = 10.0, A = 2.0, INCON X0 = 1.0.
   • TABLE for function tables.
3. Control statements: control the run and output.
   • TIMER FINTIM = 20.0, OUTDEL = 1.0, PRDEL = 0.5 set finish time, output interval and print interval.
   • PRINT Y, DY, PRTPLT, LABEL, METHOD RKS, END, STOP.

Example (decay $dy/dt=-ky$)

Y = INTGRL(Y0, DY)
DY = -K*Y
CONSTANT Y0 = 100.0
PARAMETER K = 0.3
TIMER FINTIM = 10.0, OUTDEL = 1.0
PRINT Y
END
STOP

The statement is true. In a telephone system simulation a call is a temporary entity that passes through a series of states, and a call is lost whenever it reaches the point where no free path to the called party exists.

States of a call

1. Idle / on-hook: The subscriber is not using the line.
2. Call origination (off-hook): The caller lifts the handset and the call arrives in the system. The arrival time is generated randomly (usually exponential inter-arrival times).
3. Dialling / connecting: The exchange tries to set up a path from the caller to the called party through the trunk lines (link).
4. Availability test: The simulation checks two conditions:
   • Is a link/trunk free?
   • Is the called line free (not already engaged)?
5. Connected / conversation: If both are free, the link and line are seized, the call proceeds for the holding time (a random service time) and then the resources are released. This is a completed call.
6. Busy / blocked: If the link is not available, or the called line is busy, the call cannot be set up. There is no queue in a loss system, so the caller hears a busy tone.
7. Call lost / abandoned: The blocked call leaves the system and the lost-call counter is increased. The caller may retry later as a new call.

Verification

Flow of the simulation: arrival, then test (link free AND line free). If the answer is yes, then connected, then conversation, then release, then completed. If the answer is no, the call goes to busy and is lost.

Since a call that finds either the link unavailable or the line busy has no other state to move to except the lost state, the call gets lost exactly under the stated conditions: "when the link is not available or the line is busy". Calls are lost only in this branch, and the loss probability (number lost divided by number of attempts) is a measure of performance used to decide how many links are needed.

In a discrete-event simulation the system state changes only at certain instants (events such as an arrival or the end of a service), and nothing changes between them. The next-event time advance mechanism is the way the simulation clock is moved forward.

Purpose

• To advance the simulation clock directly from the time of one event to the time of the next scheduled event, skipping the idle periods in between, instead of moving in small fixed steps.
• To save computing time, since no processing is done when nothing happens.
• To keep the model exact, because each event is executed at its precise time and not rounded to a time step.

How it works

1. Initialise the clock to 0 and place the first events in the event list (future event list), ordered by time.
2. Remove the imminent (earliest) event from the list.
3. Set the clock to the time of that event.
4. Execute the event: update the system state and statistics, and schedule any new future events (for example, an arrival schedules the next arrival and, if the server is free, a departure).
5. Repeat until the stopping condition (end time or number of events) is met.

Example: if the event list holds arrivals at $t=3$ and $t=7$ and a departure at $t=5$, the clock jumps $0 \to 3 \to 5 \to 7$.

It contrasts with the fixed-increment method, which advances by $\Delta t$ and may waste time or lose precision. Next-event advance is used in languages such as GPSS and SIMSCRIPT.

(a) Kolmogorov-Smirnov (KS) test for uniformity

The KS test checks whether a set of random numbers follows the uniform distribution on $(0,1)$ by comparing the empirical cumulative distribution $S_N(x)$ of the sample with the theoretical uniform CDF $F(x)=x$. The statistic is the largest absolute difference between them.

Hypotheses: $H_0$: the numbers are uniformly distributed; $H_1$: they are not.

Steps:

1. Rank the $N$ numbers in ascending order $R_{(1)}\le R_{(2)}\le\dots\le R_{(N)}$.
2. Compute $D^+=\max_i\left(\frac{i}{N}-R_{(i)}\right),\qquad D^-=\max_i\left(R_{(i)}-\frac{i-1}{N}\right)$
3. $D=\max(D^+,D^-)$.
4. Compare with the critical value $D_\alpha$ at the significance level. If $D \le D_\alpha$, do not reject $H_0$; otherwise reject it.

(b) Application

Sorted data, $N=6$: 0.11, 0.45, 0.54, 0.68, 0.73, 0.98.

$D^+ = 0.1033,\qquad D^- = 0.2833$ $D=\max(0.1033,\,0.2833)=0.2833$

Critical value (given) at $\alpha=0.05$: $D_\alpha = 0.565$.

Since $D = 0.2833 < 0.565$, the null hypothesis is not rejected. There is no evidence at the 5% level that the numbers differ from a uniform distribution, so the set can be accepted as uniform.

SIMSCRIPT is a general purpose, process and event oriented simulation language. Its execution is controlled by a timing routine that works on the event set (pending events) in the following cycle.

Elements

• Event notice: A record of a scheduled event holding the event type, the time at which it will occur and any attributes.
• Event set: The list of event notices ordered by increasing time of occurrence (the future event list).
• Timing routine: The built-in routine that selects and activates events.
• Event routine: The user-written code that updates the system state when the event occurs.

Execution cycle

1. Initialisation: The main program creates the entities, sets the initial state, and schedules the first events (for example the first arrival) with schedule ... now/in .. units.
2. Selection: The timing routine looks at the event set and removes the event notice with the earliest time.
3. Clock update: The simulation clock time.v is advanced to the time of this event notice.
4. Execution: The corresponding event routine is called. It changes the state variables, collects statistics, and may schedule new events (adding notices to the event set), cancel events or create/destroy entities.
5. Return: When the routine ends, control returns to the timing routine.
6. Termination test: The cycle repeats from step 2 until the event set is empty, the stop time is reached, or a stop statement is executed. Then the report is printed.

In process-oriented SIMSCRIPT, a process activates, works or waits for time, requests and relinquishes resources; the process notices are stored in the event set, and the same cycle resumes each process at the correct time.

Simple form: schedule, take the earliest event, advance the clock, execute the event routine, and repeat.

Model description

• Workers arrive every $7\pm2$ minutes (uniform between 5 and 9).
• Two clerks serve requisitions, taking $8\pm4$ minutes (uniform 4 to 12); modelled as a storage of capacity 2.
• One storekeeper, with a single facility, fills requisitions in $4\pm2$ minutes (uniform 2 to 6).
• The total time in the system for 800 requisitions is tabulated.

Block diagram (flow of blocks)

GENERATE 7,2        arrival of a worker/requisition
   |
QUEUE CLQ           join the clerk queue
   |
ENTER CLERKS        take one of the 2 clerks
   |
DEPART CLQ          leave the queue
   |
ADVANCE 8,4         clerk processes the requisition
   |
LEAVE CLERKS        free the clerk
   |
QUEUE SKQ           join the storekeeper queue
   |
SEIZE KEEPER        take the storekeeper
   |
DEPART SKQ          leave the queue
   |
ADVANCE 4,2         storekeeper fills requisition
   |
RELEASE KEEPER      free the storekeeper
   |
TABULATE TIMES      record total time in system
   |
TERMINATE 1         remove transaction, count 1

GPSS program

CLERKS   STORAGE   2
TIMES    TABLE     M1,10,10,20
         GENERATE  7,2
         QUEUE     CLQ
         ENTER     CLERKS
         DEPART    CLQ
         ADVANCE   8,4
         LEAVE     CLERKS
         QUEUE     SKQ
         SEIZE     KEEPER
         DEPART    SKQ
         ADVANCE   4,2
         RELEASE   KEEPER
         TABULATE  TIMES
         TERMINATE 1
         START     800
         END

Explanation

• GENERATE 7,2 creates transactions with inter-arrival time uniform on 5 to 9 minutes.
• The storage CLERKS of capacity 2 represents the two clerks, so up to two requisitions are processed together; others wait in the queue CLQ.
• The ADVANCE 8,4 block holds the transaction for 4 to 12 minutes.
• The single storekeeper is a facility (SEIZE/RELEASE); other transactions wait in SKQ.
• M1 is the transit time of the transaction (time since it entered the model), so TABULATE TIMES builds a frequency table of the total time for each requisition, with the first class limit at 10 and class width 10 over 20 classes.
• Each TERMINATE 1 decreases the termination counter, and START 800 ends the run after 800 requisitions have been filled. The output gives the table of times, queue statistics, and clerk and storekeeper utilisation.

(a) Need for analysis of simulation output

A stochastic simulation uses random inputs, so the output of a single run is only one sample of a random process and differs from run to run. Without proper analysis we may draw wrong conclusions. Analysis is needed:

• To estimate the true performance measures (mean waiting time, utilisation, throughput) from sample results.
• To measure the precision (variance, confidence interval) of those estimates.
• To decide the run length, warm-up period and number of replications required.
• To deal with the autocorrelation of successive observations in one run and with initial bias.
• To compare alternative system designs reliably.
• To validate the model against the real system.

(b) Estimation approach

The aim is to estimate a parameter $\theta$ (for example the mean waiting time) from output observations $X_1,\dots,X_n$.

1. Point estimation: The sample mean is used as the estimate $\bar X = \frac{1}{n}\sum_{i=1}^{n}X_i,\qquad S^2=\frac{1}{n-1}\sum_{i=1}^{n}(X_i-\bar X)^2$
2. Interval estimation: A confidence interval shows the reliability of the estimate: $\bar X \pm t_{\alpha/2,\,n-1}\frac{S}{\sqrt n}$ For example, a 95% interval means that in 95% of repeated experiments the interval would contain the true mean.
3. Terminating simulations (a natural end, such as a bank open from 9 to 5) use the method of independent replications: run $n$ replications with different random number streams, take one value $Y_j$ (the mean of run $j$) from each, and apply the formulas above to the $Y_j$, which are independent.
4. Steady-state (non-terminating) simulations require removing the initial transient by a warm-up deletion, then using one of:
   • Replication-deletion: several runs, delete the warm-up from each.
   • Batch means: one long run split into $k$ batches, whose means are treated as nearly independent observations.
   • Regenerative method: split the run at regeneration points.
5. Run length: Increase $n$ or the number of replications until the half-width of the interval is small enough, since it decreases roughly as $1/\sqrt n$.

The result is reported as a point estimate together with its confidence interval.

Bootstrapping is a resampling technique used in simulation to estimate the sampling distribution of a statistic (mean, variance, confidence interval) from a single observed data set, without assuming a particular theoretical distribution.

Procedure

1. Take the original sample $x_1, x_2, \dots, x_n$ of size $n$.
2. Draw a bootstrap sample of size $n$ with replacement from it, using random numbers to pick indices. Some observations appear more than once and some not at all.
3. Compute the statistic of interest $\hat\theta^*$ (e.g. the mean) for this resample.
4. Repeat steps 2-3 a large number $B$ of times (e.g. 1000) to obtain $\hat\theta^*_1, \dots, \hat\theta^*_B$.
5. Use these values to estimate the bias, standard error and confidence interval:

$\widehat{SE}=\sqrt{\frac{1}{B-1}\sum_{b=1}^{B}\left(\hat\theta^*_b-\bar{\theta}^*\right)^2}$

A $95\%$ percentile interval is taken from the $2.5$th and $97.5$th percentiles of the sorted $\hat\theta^*_b$ values.

Uses in simulation

• Estimating confidence intervals of output measures when the output distribution is unknown.
• Validating a model by comparing simulated and real data.
• Working with small samples where normal assumptions fail.

Advantages: simple, distribution-free, works for almost any statistic. Limitations: needs many computations, and it cannot fix a sample that is unrepresentative or too small.

Pseudorandom numbers are numbers generated by a deterministic algorithm that appear random and pass statistical tests of randomness (uniformity and independence), yet are completely reproducible from a starting value called the seed. They are used in simulation because truly random numbers are slow to produce and cannot be repeated.

Generation (linear congruential method)

$X_{i+1}=(aX_i+c)\bmod m,\qquad R_i=\frac{X_i}{m}$

Here $X_0$ is the seed, $a$ the multiplier, $c$ the increment and $m$ the modulus. For example, with $a=5$, $c=3$, $m=16$, $X_0=7$: $X_1=(35+3)\bmod 16=6$, $X_2=(30+3)\bmod 16=1$, $X_3=(5+3)\bmod 16=8$, giving $R_1=0.375$, $R_2=0.0625$, $R_3=0.5$.

Desirable properties

• Uniformly distributed on $(0,1)$.
• Statistically independent (no correlation between successive values).
• Long period before the sequence repeats.
• Fast, uses little memory.
• Reproducible, so experiments can be repeated and compared.

Tests: chi-square or Kolmogorov-Smirnov for uniformity; runs test and autocorrelation test for independence.

Drawback: the sequence eventually repeats, and a poor choice of $a, c, m$ gives short periods or patterns. Pseudorandom numbers are then converted by inverse transform or other methods into variates from other distributions.

The parasite-host model is a biological (continuous or discrete) system model describing how the populations of a host species and the parasite that lives on it affect each other. It is a predator-prey type interaction and a standard example of simulating a dynamic system with feedback.

Basic idea

• The parasite depends on the host for survival and reproduction, so parasite growth requires hosts.
• Parasites reduce host survival or reproduction.
• When hosts are abundant, parasites multiply; as parasites increase, hosts decline; then fewer hosts cause the parasite population to fall, allowing hosts to recover. This produces cyclic oscillations.

Mathematical form (Lotka-Volterra type)

Let $H$ be the host population and $P$ the parasite population:

$\frac{dH}{dt}=aH-bHP,\qquad \frac{dP}{dt}=cHP-dP$

where $a$ is the host growth rate, $b$ the rate at which parasites reduce hosts, $c$ the parasite growth gained from hosts, and $d$ the parasite death rate. The equilibrium is $H^*=d/c$, $P^*=a/b$.

A discrete version is the Nicholson-Bailey model, in which $H_{t+1}=\lambda H_t e^{-aP_t}$ and $P_{t+1}=cH_t\left(1-e^{-aP_t}\right)$.

Simulation: starting from initial values $H_0, P_0$, the equations are solved step by step (e.g. Euler's method with a small time step) to observe how the populations change and to test the effect of changing parameters. The model is applied in epidemiology, pest control and ecology.