Introduction to Modern Mathematical Modeling with R:
A User’s Manual to Train Mathematical Consultants
A Cambridge University Press book by
SSP Shen
R Code written by Dr. Samuel Shen, Distinguished Professor
San Diego State University, California, USA
Email: sshen@sdsu.edu
Compiled and Edited by Joaquin Stawsky
San Diego State University, August 2022
Chapter 8: Stochastic Models
A Nowhere Differentiable but Everywhere Continuous Model
#Wierstrass function
a <- 0.5
b <- 13
m <- 10000
n <- 100
w <- matrix(rep(0,m*n),ncol=n)
x <- seq(-2,2,length=m)
w[,1] <- (a^(+1))*cos(b*pi*x)
for (k in 2:n) w[,k] <- w[,k-1]+(a^(+k))*cos((b^k)*pi*x)
plot(x, w[,n]+1,type="l")
grid(nx = NULL, ny = NULL)
#Riemann's construction
m <- 10000
n <- 10000
wr <- matrix(rep(0,m*n),ncol=n)
x <- seq(-2,2,length=m)
wr[,1] <- sin(x)
for (k in 2:n){
wr[,k] <- wr[,k-1]+(sin((k^2)*x))/(k^2)
}
plot(x, wr[,n],type="l")
grid(nx = NULL, ny = NULL)
Brownian Motion
#Simple Brownian Motion Code in One Dimension N <- 1000
N <- 1000
dis <- rnorm(N, 0, 1);
dis <- cumsum(dis);
plot(dis, type = "l",
main = "Brownian Motion in One Dimension",
xlab = "Time", ylab = "Displacement")
grid(nx = NULL, ny = NULL)
#Simple Brownian Motion Code in 2-Dimension N <- 1000;
N <- 20
xdis <- rnorm(N, 0 ,1);
ydis <- rnorm(N, 0 ,1);
xdis <- cumsum(xdis);
ydis <- cumsum(ydis);
plot(xdis, ydis, type="l",
main ="Brownian Motion in 2-Dimension",
xlab="x", ylab = "y",
col="blue")
grid(nx = NULL, ny = NULL)
#Put arrows on the point locations for the order of particle motion
s <- seq(length(xdis)-1)
arrows(xdis[s], ydis[s], xdis[s+1], ydis[s+1],
col="red")
#Brownian Motion with Displacement Code in One Dimension N <- 1000
N <- 1000
dis <- rnorm(N, 0, 1);
at <- rpois(N,1)
for(i in 1:N){
if(at[i] != 0){
dis[i] <- dis[i]*at[i];
}
}
dis <- cumsum(dis)
plot(dis, type= "l",
main= "Brownian Motion in One Dimension with Poisson Arrival Process",
xlab="Time", ylab="Displacement",
cex.lab=1.2, cex.axis = 1.2)
grid(nx = NULL, ny = NULL)
#Brownian Motion With Displacement Code in Two Dimension N <- 1000
N <- 1000
xdis <- rnorm(N, 0,1);
ydis <- rnorm(N, 0,1);
xdis <- cumsum(xdis);
ydis <- cumsum(ydis);
at <- rpois(N,1)
for(i in N){
if(at[i] != 0){
xdis[i] <- xdis[i]*at[i]; ydis[i] <- ydis[i]*at[i] }
}
plot(xdis, ydis, type="l",
main = "Brownian Motion in Two Dimension with Poisson Arrival Process",
xlab="x Displacement" , ylab="y Displacement")
grid(nx = NULL, ny = NULL)
Dimension of Koch Curve
#TurtleGraphics package
#install.packages("TurtleGraphics")
library(grid)
library("TurtleGraphics")
turtle_init() 
turtle_forward(dist=30)
turtle_backward(dist=10) 
turtle_right(angle=90) 
turtle_forward(dist=10) 
turtle_left(angle=135) 
turtle_forward(dist=14) 
turtle_left(angle=90) 
turtle_forward(dist=14) 
turtle_left(angle=135) 
turtle_forward(dist=10)
#Koch snowflake
koch <- function(s=50, n=6) {
if (n <= 1)
turtle_forward(s)
else {
koch(s/3, n-1)
turtle_left(60)
koch(s/3, n-1)
turtle_right(120)
koch(s/3, n-1)
turtle_left(60)
koch(s/3, n-1)
}
}
turtle_init(600, 400, "error") 
turtle_do({
turtle_up()
turtle_left(90)
turtle_forward(250)
turtle_right(180)
turtle_down()
koch(500, 6)
})
Use R to Calculate the Fractal Dimension
# R method for estimating fractal dimensions for time series:
#install.packages("RandomFields")
#library(RandomFields)
#install.packages("fractaldim")
library(fractaldim)## Loading required package: abind#Standard normal time series
rf2 <- rnorm(1000)
fd.estimate(rf2, methods="variation",
plot.loglog = TRUE, col="blue") 
#This fd.estimate command yields a figure which shows the logN vs log epsilon line
#whose slope is the negative of the Hausdorff dimension we wish to compute.
#D <- 2.02
#Uniformly distributed random time series
rf2 <- rnorm(10000)
fd.estimate(rf2, methods="boxcount", plot.loglog = TRUE, col="blue") 
#D <- 1.89
#Brownian motion time series
set.seed(123)
N <- 1000
T <- 1
delt <- T/N
W <- rep(0,N+1)
for (i in 1:N){
W[i+1] <- W[i]+rnorm(1)*sqrt(delt)
}
plot(seq(1,N+1),W,type="l")
for (i in 1:N){
W[i+1] <- W[i]+rnorm(1)*sqrt(delt)
}
plot(seq(1,N+1),W,type="l")
fd.estimate(W, methods="boxcount",
plot.loglog = TRUE, col="blue") 
#D <- 1.51