Predicting CPI Performance using an Ensemble of Neural Networks

A photo of the difference in error between an ensemle of neural networks and single neural network

This project is dedicated to investigating the utility of an ensemble of recurrent neural networks used to predict fluctuations of the Consumer Price Index (CPI). Twenty technical analysis indicators are used as inputs to the network. An individual networks attempt to model a single scenario of the CPI's behavior. Models from each scenario were then averaged to produce a unique prediction. My results show that an average ensemble model consistently outperforms the prediction made by a single network. For an in depth explanation of this research, please click on this paper.

This paper explores the relationship between neural networks found in the brain, classical statistical physics, Ising type models, quantum potentials, and ensembles of recurrent neural networks. Also present in the report are several programs (out of many) that were used to clean the data and generate the desired predictions. Some of this same code can be accessed in the drop down below.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
Created on Fri Nov 26 16:43:48 2021

@author: nickdorogy

##Import statements 
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow.keras import datasets, layers, models
from sklearn.ensemble import RandomForestClassifier 
from keras.utils import to_categorical
from keras.models import Sequential
from keras.layers import Dense
from matplotlib import pyplot
import numpy as np 
from numpy import mean
from numpy import std
from numpy import array
from numpy import argmax
from numpy import tensordot
from numpy.linalg import norm
from itertools import product

##Initial data load-in and splitting data into X and y 
##X is the raw data and y is the standaridzed Consumer Price Index
X = np.genfromtxt("/Users/nickdorogy/Desktop/school/GitHub/CapstoneCode-/Capstone/Standardized/!ALL_COMBINED/ALL_COMBINED_STANDARDIZED_DATA.csv", delimiter=',', skip_header=1, usecols=(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21)) 
y = np.genfromtxt("/Users/nickdorogy/Desktop/school/GitHub/CapstoneCode-/Capstone/Standardized/Consumer Price Index/STANDARDIZEDConsumerPriceIndex.csv", delimiter=',', skip_header=1) 

##Splitting data into train and test 
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.20, random_state=33) #select size of test data (and therefore size of train)
y_test = np.reshape(y_test, (40,1)) #reshaping so everything is matrix of same dimensions 

##Fit model on dataset
def fit_model(X_train, y_train):
    baseline_model = models.Sequential()
    baseline_model.add(layers.Dense(26, activation='sigmoid')) #input 
    baseline_model.add(layers.Dense(1)) #output-- obviously want one output
    #choose optimizer, loss, & metric
                    loss=tf.keras.losses.MeanSquaredError(reduction="auto", name="mean_squared_error"),
                    metrics=['RootMeanSquaredError']), y_train, epochs=50,batch_size=1) #Select number of epochs & batch size  
    return baseline_model
##Make an ensemble prediction
def ensemble_predictions(members, X_test):
    #make predictions
    y_predict = [model.predict(X_test) for model in members] 
    #reshaping so they are matrices of the same dimensions 
    y_predict = [np.reshape(pred, (40,1)) for pred in y_predict]
    return result

##Evaluate accuracy of the model 
def evaluate_members(members, X_test, y_test):
    y_predict = ensemble_predictions(members, X_test)
    accuracy = np.abs(y_predict - y_test) 
    return accuracy

##Fit all models
n_members = 500 # number of models to be averaged 
members = [fit_model(X_train, y_train) for i in range(n_members)]

##Calculate & display error 
ensemble_score = evaluate_members(members, X_test, y_test)
single_score = np.abs(np.reshape(members[0].predict(X_test), (40,1))-y_test)
print("Ensemble Error:", sum(ensemble_score))
print("First Network's Error:", sum(single_score))   
##Summarize average accuracy of a single final model
print('Accuracy of single score || Mean: %.3f, Std Dev:, (%.3f) ||' % (mean(single_score), std(single_score)))
print('Accuracy of ensemble score || Mean: %.3f, Std Dev:, (%.3f) ||' % (mean(ensemble_score), std(ensemble_score)))

##Plot scores
plt.plot(np.abs(single_score), label="Single network error")
plt.plot(np.abs(ensemble_score), label="Ensemble of networks' error")
plt.plot(ensemble_score-single_score, label="Ensemble - single error")
plt.xlabel('Time (Months)')

New Measurements of OZ UMa

A photo of the period we estimated for OZ UMa

What are RR Lyrae Stars?

In the fall of 2020, I joined a research team studying RR Lyrae type variable stars. This class of star generally exhibits a periodic behavior that allows us to estimate luminosity, and therefore absolute magnitude according to the equation below:

M1 - M2 = -2.5 × log(L1L2)

where M1 is the absolute magnitude of star 1, M2 is the absolute magnitude of star 2, L1 is the luminosity of star 1, and L2 is the luminosity of star 2. If we examine the equation, we know that in order to calculate the absolute magnitude of star 1, we need to know its own lumonisty (easily found for RR Lyrae stars) and a second, reference stars' absolute magnitude and luminosity. So, by measuring the luminosity of a stellar object (a comparatively simple task) we are able to generate an estimate of the star’s absolute magnitude. We can then combine the absolute magnitude with our observed apparent magnitude to estimate the distance to the star as described by:

d = 10(m - M - A + 5) ∕ 5

where m is the apparent magnitude, M is the absolute magnitude, d is the distance to the star, and A is the extinction factor. For our research, we performed the absolute magnitude and distance calculation several times with different filters to ensure more accurate results. We explain this process in more detail here.

Thus, because of a simple periodic behavior inherent in RR Lyrae stars, we are able to produce an accurate distance measurement! This unique property has landed RR Lyrae stars the designation of standard candles-- objects that serve a distance reference for the surrounding stellar objects. Although perhaps relatively abstract, this places special significance on this class of stars. If we are able to locate an RR Lyrae star in a galaxy, galactic cloud, or any other celestial body, we can derive a distance measurement to that location and use it as a jumping point for increasingly farther distance calculations. With a network of such standard candles in place, we can begin charting a huge number of systems that currently have little or no recorded information. Of course, this process becomes much more challenging when applying theory to reality, but it still remains a particularly accurate and comparatively easy method for estimating distances.

Our Role with OZ UMa

Although the project is based out of Australia, we were able to conduct our work remotely in coordination with our Washington and Lee-based research advisor, David Sukow.

As members of the larger research group, our role was to select a star and carry out the aforementioned process for estimating distance. My team, consisting of two W&L classmates, selected OZ UMa has our target. OZ UMa is an RRab star (a sub-class RR Lyrae stars) situated in Ursa Major and has very little information available from its limited previous studies. Using publicly available telescope time through Las Cumbres Observatory Network, we tracked OZ UMa over a 14 night period to gather information on the stars period fluctuation. Our efforts have culminated in a number of new or updated measurements of OZ UMa’s stellar properties. Although currently awaiting publication to the Astronomy Theory, Observations, and Methods Journal, you can read more about OZ UMa and our research through this pdf.

Completed Projects

432 Park Ave. Resonant Frequency Simulation

GIF of 3D motion of building

432 Park Ave. is one of New York's City most recent additions to the city skyline. Completed in 2015, the estimated cost of the project sits between a staggering $1.25 and $3.12 billion. Although construction costs are partially attributable to its location on Billionaire’s Row, the building’s innovative design certainly influenced the astronomical cost. Measuring at a 15:1 height-to-width ratio, 432 Park Ave. is currently one of the most slender structures in the world placing it in a new class of buildings described as “Pencil Towers.”

Much like balancing a pencil vertically on its eraser, fashioning such a building requires an extreme amount of precision. But, for a pencil, you need not even blow on it to make it fall on its side-- so how does a building keep from collapsing when subjected to hurricane-force winds, earthquakes, the commotion of underground railroads, and its own internal reverberations? The simulation to the right (provided by Strand7) provides a good visualization of extreme movement. As it turns out, the solution is not to stop the building from moving, but rather to control it. For 432 Park Ave., this required innovative design procedures like “gap floors” that dilute the effects of wind.

This area of study is of extreme practical importance when designing structures and is the focus of our investigation into 432 Park Ave. Specifically, we study how the building is constructed to ensure that it is not forced to oscillate in its resonant frequency which can produce catastrophic side-effects.

A GIF of Tacoma Narrows

Natural phenomenon like earthquakes and hurricane force winds are a common and potentially dangerous influence on a building’s movement. Because these events generally oscillate at very low frequencies, when modelling their potential effects, we can limit our scope of study to the first ten natural frequencies of a system. In our case, if we know the first 10 natural frequencies of 432 Park Ave., we can compare them with historical data and machine generated simulations of natural phenomenon. This allows engineers to ensure that the structure will not have a natural frequency that is equal to a frequency produced by the environment. The GIF to the left is an example of a bridge (the Tacoma Narrows Bridge) that was forced to resonate in its natural frequency by strong gusts of wind and later collapsed.

Example of degrees of freedom

To model this system, we utilized a degree-of-freedom model where each floor added a single degree-of-freedom (see accompanying illustration to the right). Thus, with 96 total floors, including basement and “gap” floors, we have 96 degrees-of-freedom and, in turn, our model produces 96 coupled differential equations that describe the buildings motion. This allows us to consider the buildings support beams as springs, where the equation for the spring constant of a single support beam is given by:

k = 48 × (E × I) ∕ L3

Springs (or support beams) on a particular floor are considered to be in series. Therefore, once we know the spring constant of a single beam, we can multiply this value by the number of springs on that floor to find the total value of the spring constant to be used in our differential equation. In the above equation, E is Young's modulus for high compression concrete, I is the area moment of intertia for a rectangular cross section, and L is the length of each beam. For a brief overview of these values and how we calculated them, you can visit this document.

Casting our degree-of-freedom equations as a matrix allows us to calculate the modes and mode shapes of the the structure by solving for the eigenvectors and eigenvalues. When solving this system of equations, the eigenvalues correspond to the building's natural frequencies, while the eigenvectors describe the motion (or mode shapes) of the building for a specific eigenvalue. Each eigenvalue will have a corresponding eigenvector to describe its motion. As you might imagine, this is not a particularly fun equation to write out by hand, so we turned to MatLab to perform this complex computation.

To to do this, we utilized and expanded on previous MatLab code to compute our values. The drop down below contains the original code and the modifications we used to generate the eigenvectors and eigenvalues of our system. It also contains the script for the corresponding simulation rendered in MatLab.

%% script for running N story building sim.  Two examples are included.   
% 1. No initial conditions specified. In this case, modes are equally mixed with all c_n's set equal and phases all set to 0. 
% 2. Initial conditions specified.  Note that you could use the eigIC.m script to solve for C and Phi first. 
% Last modified; April 2020
%% Example with no Initial conditions specified.  
%specifying single values for m and k, NstoryBulidingSim generates the appropriate M and K matrices. 
%Other option is to fully specify M and K as matrix and pass those the main workhorse function NStoryBuildingSim(...) 

% %% 2. Example specifying initial conditions 
% Nfloors = 2; 
% m = 2e5; 
% k = 4e5; 
% Phi = [pi/2, pi/2]; 
% C = 0.2*[0.9772; -2.5583]; 
% [Omega, V] = NStoryBuildingSim(Nfloors, m, k, C, Phi);

%Our Values
m = 487350; %kg 
k = 55600000000; %N/m 
Nfloors = 96; 
close all; 
[Omega, V] = NStoryBuildingSim(Nfloors, m, k); 


% [Omega, V, simtimes, Xout, f, Fspec] = NStoryBuildingSim(Nfloors, m, k c, phi)
% Simulates building modes with NFLOORS.  
% The mass of each floor is M, the total support/sping constant for each floor is K.
% Each floor and set of support columns is assumed identical
function [Omega, V, simtimes, Xout, f, Fspec] = NStoryBuildingSim(Nfloors,m,k, c, phi)

% Initial conditions can be input with C and PHI (Nfloors x 1). 
% Plots are generated showing the mode shapes swaying in time, the displacement of each floor.  When the simulation finishes the FFT is
% displayed. Outputs allow further processing and plotting, if desired. By default, the total simulation time is equal to 3 periods of the
% largest natural period (smallest wo)
% ---- INPUTS ----
% Nfloors - number of floors in the building. Default = 3.
% c - mode mixing constants (Nfloors x 1).  Default is equal mode mixing:
%       ones(1,Nfloors)/Nfloors, e.g. 1/3*[1 1 1] 
% phi - initial phase angles. Default: zeros.
% ---- OUTPUTS-----
% Omega  - vector of natural frequencies (rad/s)
% V -  Modal matrix. Each column corresponds to a (normalized) mode shape
% simtimes - vector of times at which positions are computed
% Xout - X position of each floor vs. time. Data is rows.
% f - vector of frequencies at which frequency spectra FSPEC is computed.
% Fspec - frequency spectra. Data in rows.
% -------- EXAMPLE -----
% Nfloors = 4;
% c = [1, 0.5, 0.25, 0.75];
% NStoryBuildingSim(Nfloors, c);
% See also eig, diag
% Original Author: Jon Erickson, W&L, 
% Last updated: April 2020 

%% JE Mar 23 2010
clrwheel2 = cell2mat({'b'; 'r'; 'c';'g'; 'm'; 'k'});

Ncycles =3;
%% choose a plot style
plotoverlay = false;

%% choose initial condition coeffs;
if nargin<1 | isempty(Nfloors)
    Nfloors = 3;

if nargin<4 |isempty(c)
    c = ones(1,Nfloors)/Nfloors;

if nargin<5 | isempty(phi)
    phi = zeros(1,Nfloors);

%% define matrices and solve for eigvects and vals
%mass matrix
%assumes each floor has same mass
M = m*eye(Nfloors);

%floors 86 and above are heavier by factor of 1.89
M(86:end, :) = M(86:end, :)*1.89;
% M2 = m2*eye(Nfloors);

%stiffness matrix
%assumes all columns are the same

% if Nfloors < 86
K = diag(2*k*ones(1,Nfloors))+diag(-k*ones(Nfloors-1,1),1)+diag(-k*ones(Nfloors-1,1),-1);
K(Nfloors, Nfloors) = k;
% else
K(86:end, :) = K(86:end, :)*0.78;
% K2 = diag(2*k2*ones(1,Nfloors))+diag(-k2*ones(Nfloors-1,1),1)+diag(-k2*ones(Nfloors-1,1),-1);
% K2(Nfloors, Nfloors) = k2;
% K = [k1+k2 -k2 0 ;
%     -k2 k2+k3 -k3;
%     0 -k3 k3];

%dynamic matrix
 D =inv(M)*K;
% if Nfloors < 86
%     D =inv(M)*K;
% else Nfloors > 86
%     D2 = inv(M2)*K2;
% end

[Vall,LambdaAll] = eig(D); % V is eigenvectors in columns. D is diag matrix of eigenvals.
OmegaAll = sqrt(diag(LambdaAll )); %vector of natural frequencies (rad/s).

% grab the 10  (nmodes) loest frequencies of oscillation
Nmodes = 10;
[OmLow, IIom] = sort(OmegaAll, 'ascend');
 Omega = OmLow(1:Nmodes);
V = Vall(:, IIom);

%normalize eigenvectors
for j = 1:Nmodes
    V(:,j) = V(:,j)/norm(V(:,j));

%sort the modes according to frequency:
[Omega,I] = sort(Omega, 'ascend');
V = V(:,I);

whos V Omega c
%computer natural periods:
To = 2*pi./Omega;
DT = min(2*pi./Omega)/10;
%% plot eigenvectors/shapes.
modefig = figure(1);

% modefignum = modefig.Number;
% figure(modefig);
    set(modefig, 'units','normalized', 'position', [ 0.0469    0.2046    0.5917    0.6889]);

% set(gcf, 'position', fullfig());
xshift = 1;

% close all;
iter = 1;
simtimes = 0:DT:Ncycles*max(2*pi./Omega);  %JE 09 Mar 2018: correcte dtime interval
Xout = zeros(Nfloors,length(simtimes));

for t = simtimes
    %compute total response at time t
    %     legend(num2str(To), 'location', 'northeast')
    X = zeros(Nfloors,1);  %total response. We'll sum to find it
    whos X V omega c
    for j = 1:Nmodes
        X = X + c(j)*V(:,j)*cos(Omega(j)*t + phi(j));
    Xout(:,iter) = X;
    %then plot
    offset = zeros(1,Nmodes);
%     set(modefig, 'units','normalized', 'position', [ 0.6469    0.5046    0.2917    0.3889]);
    for j = 1:Nmodes+1
        if j==1
            %             offset(j) = max(V(:,j));
            %             offset(j) = max(V(:,j))+max(V(:,j-1));
            hold on;
            grid on;
        %% Plot either overlaid or separate modes
        if j<=Nmodes
            clrJ = 1+mod(j,length(clrwheel2));
            % this plot command overlays the modes as they are plotted.
            if plotoverlay
                plot([0;V(:,j)]*cos(Omega(j)*t+phi(j)), 0:Nmodes, [clrwheel2(clrJ),'o:'], 'markerfacecolor', clrwheel2(j));
                plot([0;V(:,j)]*cos(Omega(j)*t+phi(j))+(j-1)*xshift, 0:Nfloors, [clrwheel2(clrJ),'o:'], 'markerfacecolor', clrwheel2(clrJ));
                %                 plot([0;V(:,j)]*cos(Omega(j)*t)+xshift*(j-1), 0:Nmodes, [clrwheel2(j),'o:'], 'markerfacecolor', clrwheel2(j));
            if plotoverlay, plot([0;X], 0:Nfloors,'ko-', 'markerfacecolor', 'k'); else
                plot([0;X]+xshift*(j-1), 0:Nfloors,'ko-', 'markerfacecolor', 'k');end
    xlabelcell = {};
    for nn = 1:Nfloors
        xlabelcell{nn} = num2str(nn);
    xlabelcell{Nfloors+1} = 'Total';
    set(gca, 'xticklabel', xlabelcell, 'xtick', xshift*[0:Nmodes]);
    %     set(gca, 'xticklabel', {'1'; '2'; '3'; 'Total'}, 'xtick', xshift*[0:Nmodes]);
    if ~plotoverlay
        xlim(1.05*[-xshift, xshift*(Nmodes+1)]);
        xlim(1.05*xshift*[-1 1]);
    ylim([0 Nmodes+1]);
    title(sprintf('Time = %2.3f s', t));
    if iter==1; pause; end
    %% plot the position of each mass over time.
    set(gcf, 'units', 'normalized', 'position', [ 0.6443    0.03    0.2917    0.85]);
    for kk = 1:Nmodes %top floor is actually last mode.
        if kk ==1, xlabel('Time (s)'); end
        subplot(Nmodes,1, Nmodes-kk+1); hold on; xlim([0 simtimes(end)]);
        plot(simtimes(iter), Xout(kk,iter), 'ko-', 'markersize', 4);
        ylabel(sprintf('x_%d', kk));
        grid on;
        xlabel('Time (s)')
    iter = iter+1;
wo = 2*pi./To;
legend(num2str(To), 'location', 'northeast')

%% plot the FFT for each floor
L = length(Xout);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(Xout',NFFT)/L;
Fs = 1/DT;
f = Fs/2*linspace(0,1,NFFT/2+1);

set(gcf, 'units', 'normalized', 'position', [0.3432    0.0417    0.2917    0.3889]);
% Plot single-sided amplitude spectrum.

for ii = 1:Nmodes
    Fspec(:,ii) = 2*abs(Y(1:NFFT/2+1,ii));
    subplot(Nmodes,1,Nmodes-ii+1); %plot top floor on top. bottom floor on bottom
    axis tight
    % Fspec = abs(Y);
    title(sprintf('FFT of x%d(t)',ii))
    xlabel('Frequency (Hz)')

This process can only be described as a first order approximation at its very best. To keep this project from getting out of hand, we do not consider particular construction designs (namely outriggers) or the effects of tuned mass dampers (TMDs). Moreover, we approximate construction design and materials to the best of our knowledge. To read more about 432 Park Ave. and degree-of-freedom systems visit the links below:

Ballistic Motion Calculator

I created this ballistic motion calculator to serve as the final project for the Washington and Lee Physics Department Modelling and Simulations in Physics course. The intended (and successful!) outcome of the project was to merge my physics and coding knowledge to develop a program capable of calculating and displaying the trajectory of a projectile under user-provided conditions.

To provide a sense for the breadth of the project, listed below are some of the conditions the calculator is able to account for:

Stillshot of 3D ballistic motion animation
  1. Elevation
  2. Temperature
  3. Exit velocity/muzzle velocity
  4. Air density
  5. Projectile mass
  6. Projectile cross area
  7. Projectile diameter
  8. Projectile ballistic coefficient
  9. Target distance
  10. Wind speed and direction

These inputs were then supplied to my python program that performed the appropriate calculations to describe the projectiles motion. Drawing upon the Euler-Lagrange equations, my program is designed to update the trajectory of the projectile in discrete, small time steps. The picture below is a stillshot of the animation generated by the program as it follows the projectiles path in times steps of dt = 0.05. At the sacrifice of computation time, we can use increasingly smaller time steps to produce increasingly more accurate outputs. However, in this case, A time step of 0.05 is sufficient.

My first challenge was to acquire the relevant inputs. To do this, I used the tkinter GUI interface to prompt the user for original inputs. However, I was careful to design the input request to ensure that potential future updates would allow for these inputs to be provided by external sources. For example, a modified Raspberry Pi would capable of providing the majority of these inputs for real-time, accurate data.

Next, I had to determine the appropiate equations to account for variables such as drag, wind, elevation, etc. Because of the difficulty associated with accurately accounting for drag, I opted to request the user for the projectile's ballistic coefficient. Generally, this value is known or even listed by the manufacturer and is much better suited for performing calculations on a wide variety of projectiles that can behave quite uniquely. To consider the effects of wind, I merely needed to decompose the vector input into constiuent components (x,y, and z direction) and include these components in the equation to update velocity. This means that our equations for velocity will use the relative velocity of the projectile.

Once the variable equations are constructed, the most important step can be included; the main equations to update motion! As previously mentioned, I draw upon the Euler-Lagrange equations to create a system of second order differential equations that will update the acceleration, and thus the velocity and position, at a time step of 0.05 (units depend on the users' inputs).

Finally, we can animate the motion of our projectile! To do this, I utilize the MatPlotLib Python package. This allowed me to create a 3-D environment that tracks and plots the motion of the projectile with the appropiate forces influencing its motion.

Washington and Lee Network Research

A social network is a group of people, each of whom is connected with a subset of the others. Such a network is represented using vertices (to represent each person) that are joined by edges (representing some sort of connection).

These social networks have been a topic of interest in the social sciences for over fifty years. Models have been developed to predict the spread of disease through a group of people, the relative connectivity of random individuals, and a variety of other important aspects of social networks.

Our project was aimed at building a social network based on scientific collaboration at an elite liberal arts college. Specifically, we focused on published researchers in the sciences at Washington and Lee University. Two of these researchers are considered connected if they have ever authored a paper together. This parameter is derived from the canonical paper from The Santa Fe Institute, "The Structure of Scientific Collaboration Networks." We restricted our data to Washington and Lee’s “hard sciences” departments: biology, chemistry, geology, physics, and biochemistry. This is in response to the tendency of social science professors to abstain from collaboration.

We subsequently compare these findings to Newman’s conclusions, a well-known application of the theory regarding the average degree of connectivity between individuals. Our result establishes a connection between a small, elite liberal arts college and the research tendencies of its faculty members. Future research could be aimed at extending this research to other similar-sized institutions for comparison and analysis.

For a summary of our results, please follow this link.

COVID-19 Model


In late February of 2020, I read an academic paper attempting to model the dissemination of a new virus outbreak in Wuhan, China. Because we had just finished a unit on Monte Carlo Markov Chain modelling in my intermediate physics course, I forwarded this paper to my professor and suggested creating an SIR model to compare against the paper’s results. After discussing the article, our class split into several teams and created a python-based program that could effectively simulate the spread using an SIR model similar to the method outlined in the paper.

After discussing the article, the class split into teams to create a python-based program that could simulate the spread of the disease using an SIR model similar to the one we'd read about. Attached below is the python script for our model, along with the short write-up we completed during our class period.

Python Script

The collapsible box below contains a portion of the code we created to develop our COVID-19 SIR model.

    Created in March 2019
    Written in Python 3

    @author:  Anthony Lorson, Elise Baker, Andrew Sontchi, Nicholas Dorogy (based on a code from Oregon State)
    # This loads some pacakges that have arrays etc and the ODE integrators
    import scipy, scipy.integrate
    import pylab

    # Initial conditions
    R0 = 7400000 # number of susceptible people 
    C0 = 53 # number of latent people 
    L0 = 53 # number of infectious people 
    I0 = 2 # number of people removed (either overcome the infection, or died from it)
    alpha = .80 # infectious period
    beta = .60 # incubation period
    gamma = 0 # number of people entering 
    delta = 50000 # number of people leaving 

    Y0 = [ R0, C0, L0, I0 ]

    #ask user how long they would like to run the trial
    tMax = input("Please enter how many time steps you would like the simulation to run(in integer form);")
    tMax = int(tMax)
    # Time vector for solution (divide steps into 1,001 equal intervals)
    T = scipy.linspace(0, tMax, 1001)

    # This defines a function that is the right-hand side of the ODEs
    def rhs(Y, t, alpha, beta, gamma):
    SIR model.

    This function gives the right-hand sides of the ODEs.

    # Convert vector to meaningful component vectors
    R = Y[0]
    C = Y[1]
    L = Y[2]
    I = Y[3]

    N = R + C + L + I

    # The right-hand sides
    dR = -R/N * (2/alpha)*L + gamma - (gamma/N)*R
    dC = R/N * (2/alpha)*L - C/beta - (gamma/N)*C
    dL = C/beta - L/alpha - (gamma/N)*L
    dI = L/alpha + (gamma/N)*(L + R + C)

    # Convert meaningful component vectors into a single vector
    dY = [ dR, dC, dL, dI]

    return dY

    # Integrate the ODE
    solution = scipy.integrate.odeint(rhs,
                                args = (alpha, beta, gamma))
    R = solution[:, 0]
    C = solution[:, 1]
    L = solution[:, 2]
    I = solution[:, 3]

    N = R + C + L + I

    # Make plots

    # Load a plotting package

    pylab.plot(T, R / N, 'r',
        T, C / N, 'purple',
        T, L / N, 'b', 
        T, I / N, 'g')


    pylab.legend([ 'Suceptible', 'Latent', 'Infectious', 'Removed'])

    # Actually display the plot


Attached here is the summary of our work that we wrote to accompany the model's results and to provide further description of our work. As you can see, this model does not forecast global dissemination of the disease, but instead graphs the number of susceptible, infected, and recovered individuals as a function of time for Hong Kong specifically. Our data was limited by the availability of travel information for the city, and thus, is based on the previous year's reported statistics.

When first creating this model, we were confused by the extremely swift growth rate over such a short period of time. So, we amended our model to use different initial variables that gave a somewhat more satisfactory graph. As it turns out, we should have paid more attention to our original results!