Deltat<- 60 	# minutes
lambda1<- 1.
set.seed(1)    #  reset random numbers
# vector of the number of accepted arrivals:
ev<-  0
# vector of the number of arrivals if it were a homogeneous process:
evt<- 0
# vector of values lamda(t)/lambda1:
evs<- 0
t<-   0
while (t<Deltat) {      # arrival times generation
t<- t + rexp(1,lambda1)
evt<- c(evt,t)   # arrival times for a homogeneous process
fl=sin(t)/lambda1   # lambda(t)/lambda1
evs<- c(evs,fl)  # fl values, but not all are accepted
if (runif(1) < fl) ev<- c(ev,t)  #  <-----  crucial point
}       # ending loop while
length(ev)
length(evt)
plot(evt,0:(length(evt)-1),type="s",ylim=c(-4,length(evt)),xlab="time (min)",
ylab="No. arrivals",lty=3,col="black",cex.lab=1.1,font.lab=3,lwd=2)
curve(sin(x),0,Deltat,add=T,lty=2,lwd=1.6)
lines(ev,0:(length(ev)-1),type="s",lwd=2.2,col="black")
## Code_3_5.R
# Nonhomogeneous Poisson process is a process with intensity lambda not constant,
# but depending on time: lambda(t).
# It is a not-stationary increments process,
# not depending only on the length interval, but also on its time position.
# The simulation is performed by sampling from a homogeneous process,
# arrivals occur according the time function. Here: lambda(t) = sin(t)
# lambda1 (homogeneous process) is known and lambda(t) <= lambda1, all t
Deltat<- 60 	# minutes
lambda1<- 1.
set.seed(1)    #  reset random numbers
# vector of the number of accepted arrivals:
ev<-  0
# vector of the number of arrivals if it were a homogeneous process:
evt<- 0
# vector of values lamda(t)/lambda1:
evs<- 0
t<-   0
while (t<Deltat) {      # arrival times generation
t<- t + rexp(1,lambda1)
evt<- c(evt,t)   # arrival times for a homogeneous process
fl=sin(t)/lambda1   # lambda(t)/lambda1
evs<- c(evs,fl)  # fl values, but not all are accepted
if (runif(1) < fl) ev<- c(ev,t)  #  <-----  crucial point
}       # ending loop while
length(ev)
length(evt)
plot(evt,0:(length(evt)-1),type="s",ylim=c(-4,length(evt)),xlab="time (min)",
ylab="No. arrivals",lty=3,col="black",cex.lab=1.1,font.lab=3,lwd=2)
curve(sin(x),0,Deltat,add=T,lty=2,lwd=1.6)
lines(ev,0:(length(ev)-1),type="s",lwd=2.2,col="black")
## Code_3_3.R
# mean recurrence time
# stationary distribution = (0.4933,0.2800,0.2267)
# Weather in the Land of Oz
# States: rain (0), sun (1), snow (2)
markov<- function(x0,n,x,P) {   # starting function
s<- numeric()
s[1]<- x0              # initial state
raw<- which(x==x0)
for(i in 2:n) {
s[i]<- sample(x,1,P[raw,],replace=T)
raw<- which(x==s[i])
}
return(s)
}  # ending function
set.seed(5)         #  reset random numbers
n<- 200             # number of steps
# transition matrix:
P<-matrix(c(0.6,0.1,0.3,0.3,0.5,0.2,0.5,0.4,0.1),3,3,byrow=T)
x<-c(1,2,3)
x0<- 1
s<-markov(x0,n,x,P)
t<-c(0:(n-1))
r<-which(s==1)
d<-diff(r)  # difference between successive values
# for instance: r = 3,7,12,13, diff(r) = 4,5,1
m<-mean(d)
m                # theory 1/0.4933 = 2.03
r<-which(s==2)
d<-diff(r)
m<-mean(d)
m                # theory 1/0.2800 = 3.57
r<-which(s==3)
d<-diff(r)
m<-mean(d)
m                # theory 1/0.2267 = 4.41
## Code_3_4.R
#  Birth-death process simulation
set.seed(2)        # reset random numbers
nhistories<-  4    # number of simulated histories
Deltat<- 3         # the process is followed no more than Deltat
lambda<- 0.5    # rate of birth
mu<-     0.9    # rate of death
n.init <- 20    # initial population size
time<-matrix(,nhistories,200)
population<-matrix(,nhistories,200)
ls<- numeric()
lp<- numeric()
###  2 loops: most external loop on histories,
#             inner loop generates births/deaths
for(l in 1:nhistories)      {        # starting loop on histories
ev<- 0
np<- n.init
n<- np
t<- 0
while(t<Deltat & np>0) {    # population may become extinct
t<- t + rexp(1,(lambda+mu)*n)
# if 1: birth, if 0: death
bin<-rbinom(1,1,(lambda/(lambda+mu)))
np<-ifelse (bin==1, np<- np+1, np<- np-1)
ev<- c(ev,t)
n <- c(n,np)
}             # ending while loop
ls[l]<-length(ev)
time[l,1:ls[l]] <- ev
lp[l]<-length(n)
population[l,1:ls[l]] <-n
}          # ending loop on histories
## # uncomment the following lines to print the transition instant
# and the number of individuals over time for each history
# time[1,1:ls[1]]
# time[2,1:ls[2]]
# time[3,1:ls[3]]
# time[nhistories,1:ls[nhistories]]
# population[1,1:ls[1]]
# population[2,1:ls[2]]
# population[3,1:ls[3]]
# population[nhistories,1:ls[nhistories]]
par(mfrow=c(2,2))
plot(time[1,1:ls[1]],population[1,1:ls[1]],type="s"
,xlab="time",ylab="population")
abline(v=Deltat,lty=3)
plot(time[2,1:ls[2]],population[2,1:ls[2]],type="s"
,xlab="time",ylab="population")
abline(v=Deltat,lty=3)
plot(time[3,1:ls[3]],population[3,1:ls[3]],type="s",
xlab="time",ylab="population")
abline(v=Deltat,lty=3)
plot(time[nhistories,1:ls[nhistories]],
population[nhistories,1:ls[nhistories]],type="s",
xlab="time",ylab="population")
abline(v=Deltat,lty=3)
## Code_6_1.R
#White noise
# for reproducibility: change seed to change noise
set.seed(3)
# create noise time series
w <- rnorm(1000)
plot(as.ts(w),ylab="White noise")
# verify w is zero mean
mean(w)
var(w)
hist(w,main="",col="white")
# verify w is uncorrelated
acf(w,main="",lag=20,ylim=c(-0.2,1))
# verify w is really white
w.spec <- spectrum(w,log="no",span=200,plot=FALSE)
plot(w.spec$freq,w.spec$spec,t="l",ylim = c(0,1.2),
xlab="frequency (cycles/sample interval)",ylab="spectral density")
## Code_6_2.R
# MA(1) positive coefficient
# for reproducibiity
set.seed(13)
# generate white (Gaussian) noise
w <- rnorm(100)
x <- rep(0,100)
# generate MA(1)
x[1] <- w[1]
for (t in 2:100) x[t] <- w[t] + 0.5*w[t-1]
# plot time series
plot(as.ts(x),ylab="MA(1) process")
# plot autocorrelation function
acf(x,main="")
# plot spectral density
x.spec <- spectrum(x,log="no",span=30,plot=FALSE)
plot(x.spec$freq,x.spec$spec,t="l",ylim = c(0,3),
xlab="frequency (cycles/sample interval)",ylab="spectral density")
# MA(1) negative coefficient
# for reproducibiity
set.seed(13)
# generate white (Gaussian) noise
w <- rnorm(100)
x <- rep(0,100)
# generate MA(1)
x[1] <- w[1]
for (t in 2:100) x[t] <- w[t] - 0.5*w[t-1]
# plot time series
plot(as.ts(x),ylab="MA(1) process")
# plot autocorrelation function
acf(x,main="")
# plot spectral density
x.spec <- spectrum(x,log="no",span=30,plot=FALSE)
plot(x.spec$freq,x.spec$spec,t="l",ylim = c(0,2),
xlab="frequency (cycles/sample interval)",ylab="spectral density")
## Code_6_2.R
# MA(1) positive coefficient
# for reproducibiity
set.seed(13)
# generate white (Gaussian) noise
w <- rnorm(100)
x <- rep(0,100)
# generate MA(1)
x[1] <- w[1]
for (t in 2:100) x[t] <- w[t] + 0.5*w[t-1]
# plot time series
plot(as.ts(x),ylab="MA(1) process")
# plot autocorrelation function
acf(x,main="")
# plot spectral density
x.spec <- spectrum(x,log="no",span=30,plot=FALSE)
plot(x.spec$freq,x.spec$spec,t="l",ylim = c(0,3),
xlab="frequency (cycles/sample interval)",ylab="spectral density")
# MA(1) negative coefficient
# for reproducibiity
set.seed(13)
# generate white (Gaussian) noise
w <- rnorm(100)
x <- rep(0,100)
# generate MA(1)
x[1] <- w[1]
for (t in 2:100) x[t] <- w[t] - 0.5*w[t-1]
# plot time series
plot(as.ts(x),ylab="MA(1) process")
# plot autocorrelation function
acf(x,main="")
# plot spectral density
x.spec <- spectrum(x,log="no",span=30,plot=FALSE)
plot(x.spec$freq,x.spec$spec,t="l",ylim = c(0,2),
xlab="frequency (cycles/sample interval)",ylab="spectral density")
## Code_6_2.R
# Moving-average process
# MA(1) positive coefficient
# for reproducibiity
set.seed(13)
# generate white (Gaussian) noise
w <- rnorm(100)
x <- rep(0,100)
# generate MA(1)
x[1] <- w[1]
for (t in 2:100) x[t] <- w[t] + 0.5*w[t-1]
# plot time series
plot(as.ts(x),ylab="MA(1) process")
# plot autocorrelation function
acf(x,main="")
# plot spectral density
x.spec <- spectrum(x,log="no",span=30,plot=FALSE)
plot(x.spec$freq,x.spec$spec,t="l",ylim = c(0,3),
xlab="frequency (cycles/sample interval)",ylab="spectral density")
