assign5.R

nesau — Mar 6, 2014, 11:25 AM

#### Vars ####
delta= 0.075
u = 0.02

#### Functions ####
# exp^(-delta*n)
v <- function(n,moment=1) {
  return(exp(-delta*n*moment))
}

# P(T(x)=t)
fxt <- function(t) {
  return(u*exp(-u*t))
}

tpx <- function(t) {
  return(exp(-u*t))
}

tqx <- function(t) {
  return(1-tpx(t))
}

udeferredtqx <- function(u,t) {
  return(tpx(u) - tpx(u+t))
}


derivative <- function(f,x) {
  epsilon = 0.00000001
  deriv = (f(x+epsilon)-f(x))/epsilon
  return(deriv)
}

# P(Z<z)
cdfZ <- function(z) {  
  cdf = 0

  if(z>=0) 
    cdf = ifelse (z<3000*v(20), cdf+tpx(20), cdf+tpx(20))

  if (3000*v(20)<=z) 
    cdf = ifelse (z< 2000*v(5), cdf+udeferredtqx(-log(z/3000)/delta, (20 + log(z/3000)/delta)), 
              cdf + udeferredtqx(10.41,9.59))

  if (2000*v(5) <= z)
    cdf = ifelse (z < 2000, cdf+udeferredtqx(-log(z/2000)/delta,(5+log(z/2000)/delta)) + 
                    udeferredtqx(-log(z/3000)/delta,(10.41+log(z/3000)/delta)), 
                    cdf + tqx(5) + udeferredtqx(5.41,5) )

  if (2000<= z) 
    cdf = ifelse (z<= 3000*v(5), cdf + udeferredtqx(-log(z/3000)/delta,(5.41+log(z/3000)/delta)),
                  cdf + udeferredtqx(5,0.41) )

  return(cdf)
}

pdfZ <- function(z) {
  if(z==0) return (tpx(20))
  return(derivative(cdfZ,z))
}

nums = seq(-10,2200,length=1000)
cdfnums <- numeric(0)
pdfnums <- numeric(0)

for (num in nums) {
  cdfnums = c(cdfnums,cdfZ(num))
  pdfnums = c(pdfnums, pdfZ(num))

}

plot(nums, ylim=c(0,1), cdfnums , ylab = "P(Z<=z)", xlab ="z", main="Cdf of Z", col="light blue")

# note that there is a proability mass of 20px for the next plot at z=0 which doesnt show up
plot(nums, pdfnums, ylab= "P(Z=z)", xlab = "z", main = "Pdf of Z", col="light green")