ST 321 Discrete Event
Simulation Approach
This lab prepares you for the next homework assignment (homework 7), which is attached at the end of this handout. The main emphasis of this lab session is on simulating stochastic models using the discrete event simulation approach (Chapter 6 of Ross).
Task 1: In class we have seen that, to generate n exponential random variables with rate lambda, we can use the following R statement: -1/lambda*log(runif(n)). Use the following codes to compare this with the R command rexp(n,lambda).
n=1000
lambda=5.6
r.values=rexp(n,lambda)
my.values=-1/lambda*log(runif(n))
par(mfrow=c(2,1)) # to split the graphic window into two
hist(r.values,probability=T,nclass=20,ylim=c(0,6),xlim=c(0,1.5))
hist(my.values,probability=T,nclass=20,ylim=c(0,6),xlim=c(0,1.5))
Task
2: Enter the following two R functions for our inventory example (Ch 6.5 of
Ross):
inventory<-function(s,S,x0,T) {
t<-0
x<-x0
y<-0
C<-0
H<-0
R<-0
infinity<-10000000000000
h<-0.1
r<-3
L<-2
lambda<-2
c<-function(y){return(y^(0.9)*1)}
t0<-rexp(1,lambda)
t1<-infinity
while (t<T) {
if (t0<t1) {
H<-H+(t0-t)*x*h
t<-t0
D<-rgeom(1,0.7)+1
w<-min(D,x)
R<-R+w*r
x<-x-w
if
((x<s)&(y==0)) {
y<-S-x
t1=t+L
}
t0=t+rexp(1,lambda)
#cat(t,x,y,R,C,H,(R-C-H)/T,"\n")
} else {
H<-H+(t1-t)*x*h
t<-t1
C<-C+c(y)
x<-x+y
y<-0
t1<-infinity
#cat(t,x,y,R,C,H,(R-C-H)/T,"\n")
}
}
profit<-(R-C-H)/T
return(profit)
}
estimateprofit<-function(s,S,x0,T,n=100)
{
estprofit<-0
for (i in 1:n) {
temp<-inventory(s,S,x0,T)
estprofit<-estprofit+temp
}
return(estprofit/n)
}
Use the above functions to estimate the profit for s=10,S=20,x0=15 and T=365:
estimateprofit(10,20,15,365)
Repeat with (s=5,S=20) and (s=15,S=20). Out of the three combinations of (s,S), which one gives you the highest estimated profit?
Homework 7 (Due
date:
This homework assignment is mainly concerned with R programming. You should do this homework assignment in a group of two. You will need to provide a hardcopy of your answers (program listings, numerical answers etc), as well as emailing your programs in Question 4 and the Bonus Question to us. This homework carries a heavier weight than other homeworks.
Question 1: (10 points) Using the inverse transform method, write an R statement to generate 3 random variables having probability density function f(x)=x^2/9,0<x<3. Hint: read Example 5a of Ross.
Question 2: (10 points) For the inventory example: with the same set-up as above, try your best to find the best combination of (s,S).
Question
3: (20 points) Suppose now that the cost c(y) of y unit is c(y)=2y+1, and that
the inventory holding cost per unit h is h=0.2. Modify the R function inventory and estimate
the profit when s=10,S=20,x0=15 and T=365.
Question
4: (60 points) Write an R function to implement the
insurance company model discussed in class (see also Ch 6.6 of Ross). For the claim distribution F, use Poisson with mean 5 (i.e., when
generate Y, use Y<-rpois(1,5)). Run your program a large number of times to
estimate the probability that the company will survive for the following
settings: a0=300,n0=10,c=7,T=10,nu=50,mu=0.5 and lambda=1.5 (you should
get an estimated probability around 0.20).
Email your R function to our grader Hui Zhang
at zhang@stat.colostate.edu. Your
program will be further tested by our grader.
Bonus
Question: (10 points) Write an R function to simulate
a homogeneous Poisson process with rate lambda. This time you should email your
function to Thomas Lee at