Mathematical Modeling of Complex Systems

Dr. Courtney Brown

Assignment #1

In this assignment, you will have fun with graphics while you investigate the various potential behaviors of the first-order linear difference equation with constant coefficients, as per our class discussion. This will get you started working with R, and you see that you can easily produce high quality "print ready" graphics. You are going to produce six pairs of plots and a regression analysis.

Use the program below by modifying two of the parameters (the slope and the intercept) as well as the initial condition to produce the pairs of plots for each difference equation. Within each pair, you will produce a trajectory plot, as well as the plot of the first differences. All of this is done in the R code below. The equation that you will be working with is Y(t+1) = aY(t) + b, where parameter "a" is the slope and parameter "b" is the intercept of the difference equation. The trajectory plot has Y(t) on the vertical axis and time (t) on the horizontal axis. Notice that this trajectory is (usually) not a line. Then you need to plot the first differences for each equation, which is Y(t+1) onto Y(t). You then need to run a regression of Y(t+1) onto Y(t). Take a look at the slope and intercept that comes back from the regression, and then look again at the plot of the first differences. Think about this.

Why do you need six pairs of plots? You need to produce difference equation trajectories that have the following behaviors: (1) unbounded increasing monotonic, (2) unbounded decreasing monotonic, (3) bounded increasing monotonic (convergent), (4) bounded decreasing monotonic, (5) convergent bounded oscillatory, (6) divergent unbounded oscillatory, (7) oscillatory and finite. Play with the slope, intercept, and initial condition in the code below to obtain the various behaviors. Play with different values to see what happens.

Describe what you have learned. Five typed pages plus figures should do it. Hand it all in.

NOTE: You may be interested to know that the concept of this assignment was originally developed by Professor John Sprague.

GRADUATE STUDENTS ONLY: Try adding a forced oscillator to the first-order difference equation and compare the results obtained without the forced oscillator. You may want to use m[sin(pt)], where m controls the magnitude of the sinusoidal oscillator and p controls the period. This will act as an imput to the first-order system. What have you learned from this? Include your thoughts in your write-up.

Here is the R code that you will use to get started. Just cut and paste it into your own R program file. Run it, and then modify it to give you the desired results.

y1 <- .3 # This is the initial condition
y2 <- 0
t <- 0
a <- -.6 # This is the slope
b <- .3 # This is the intercept
ystar <- b/(1-a) # This is the equilibrium value
timeserieslength <- 10
for (i in 1:timeserieslength) {
y2[i] <- (a*y1[i])+b
t[i] <- i
if (i < timeserieslength) y1[i+1]=y2[i]
}

mydata <- cbind(y1, y2, t)
parameters <- cbind(a, b, ystar)
parameters
mydata
plot(t, y1, type="o", ylab="Y", xlab="time", pch=20, lwd=2)
title(main = list("Figure 1: Y over time", cex=1.5,
col="black", font=1))

windows()

diff.model <- lm(y2 ~ y1)
summary(diff.model)
plot(y1, y2, xlab="", ylab="", xlim=c(0,.4), ylim=c(0,.4), pch=19)
title(xlab="Y", ylab="Y(t+1)", main="Figure 2: Plot of the first differences", cex=1.5,
col="black", font=2)
abline(diff.model, lwd=2)

 

For those interested in knowing how this assignment is done using SAS, the SAS code is presented below.

GOPTIONS lfactor=10 hsize=6 in vsize=6 in horigin=1 in vorigin=1 in;
options nocenter;
**********************************************************;
** CLASS, NOTE THAT IF YOU BEGIN A LINE WITH AN ASTERISK *
** THEN YOU CAN PUT NOTES IN YOUR PROGRAM FILES. THIS IS
** LIKE A COMMENT CARD IN SPSS. HOWEVER, REMEMBER THE FINAL SEMICOLON;
***********************************************************;
* NOTE THAT I INDENT SOME STATEMENTS. THIS
* IS JUST FOR NEATNESS.;
***********************************************************;
* COPYRIGHT (c) Courtney Brown 2004, All Rights Reserved;
* Permission granted to use this file and computer code for any nonprofit and
* educational purposes, including classroom instruction.
* No further permission required.
* Please cite source as "From www.courtneybrown.com";
OPTIONS NOCENTER;
DATA;
Y1=.3; * This is the initial condition;
A=-.6; * This is the slope;
B=.3; * This is the intercept;
YSTAR = B/(1-A); * This is the equilibrium value;
DO T=1 TO 10;
Y2=(A*Y1)+B;
OUTPUT;
Y1=Y2;
END; run;
SYMBOL1 I=JOIN color=black;
SYMBOL2 I=none font=simplex;
PROC PRINT;run;
PROC GPLOT;
PLOT Y2*T;
TITLE 'Figure 1a: A Time-Series Plot with a=-.6, b=.3, Y(0)=.3';
PROC REG;
MODEL Y2=Y1;
PROC GPLOT;
axis1 color=black order=0 to 1 by 0.1 minor=none
value=(h=1.5 f=swissb c=black)
label=(h=1.3 a=90 r=0 f=swissb c=black 'Y(t+1)');
axis2 color=black order=0 to 1 by 0.1 minor=none
value=(h=1.5 f=swissb c=black)
label=(h=1.3 f=swissb c=black 'Y(t)');
PLOT Y2*Y1 /
vaxis=axis1 haxis=axis2 vminor=0 hminor=0;
TITLE 'Figure 1b: A First-Difference Plot with a=-.6, b=.3, Y(0)=.3';
run;
quit;