[R] Weibull maximum likelihood estimates for censored data

From: Nadia Theron <NadiaT_at_investec.co.za>
Date: Tue, 08 Apr 2008 12:14:45 +0200


Hello!  

I have a matrix with data and a column indicating whether it is censored or not. Is there a way to apply weibull and exponential maximum likelihood estimation directly on the censored data, like in the paper: Backtesting Value-at-Risk: A Duration-Based Approach, P Chrisoffersen and D Pelletier (October 2003) page 8?  

The problem is that if I type out the code as below the likelihood ratio is greater than one.

> Interest

     D C

1 17 1

2 10 0

3 15 0

4 2 0

5 42 0

6 53 0

7 193 0

8 11 0

9 2 0

10 8 0

11 12 1  

library(stats4)

dur_ind_test = function (CDMatrix) # Matrix with durations and censores

{

lLnw <- function(b){

D = CDMatrix

NT = nrow(D)

a =((NT-D[1,2]-D[NT,2])/ sum(D[,1]^b))^(1/b)

f = sum(log((a^b)*b*(D[2:(NT-1),1]^(b-1))*exp(-((a*D[2:(NT-1),1])^b))))

fd1 = (a^b)*b*(D[1,1]^(b-1))*exp(-(a*D[1,1])^b)

fdn = (a^b)*b*(D[NT,1]^(b-1))*exp(-(a*D[NT,1])^b)

S1 = exp(-(a*D[1,1])^b)

SN = exp(-(a*D[NT,1])^b)

-(D[1,2]*log(S1)+(1-D[1,2])*log(fd1)+ f + D[NT,2]*log(SN)+ (1-D[NT,2])*log(fdn))

}

lLne <- function(A){

D = CDMatrix

NT = nrow(D)

b=1

f = sum(log(A*b*(D[2:(NT-1),1]^(b-1))*exp(-(A^(1/b)*D[2:(NT-1),1])^b)))

fd1 = A*b*(D[1,1]^(b-1))*exp(-(A^(1/b)*D[1,1])^b)

fdn = A*b*(D[NT,1]^(b-1))*exp(-(A^(1/b)*D[NT,1])^b)

S1 = exp(-(A^(1/b)*D[1,1])^b)

SN = exp(-(A^(1/b)*D[NT,1])^b)

lLw = D[1,2]*log(S1)+(1-D[1,2])*log(fd1)+ f + D[NT,2]*log(SN)+ (1-D[NT,2])*log(fdn)  

-(D[1,2]*log(S1)+(1-D[1,2])*log(fd1)+ f + D[NT,2]*log(SN)+ (1-D[NT,2])*log(fdn))

}  

fit1 <- mle(lLnw,start = list(b = 0.5)) # Estimate parameters using ml

fit2 <- mle(lLne,start = list(A = 0.02))

Lw <- lLnw(coef(fit1))                      # Maximum log likelihood :
Weibull            

Le <- lLne(coef(fit2))                      # Maximum log likelihood :
Exponential
LR0 <- (Le/Lw)                        # Likelihood ratio with duration
sample  

NSimM <- cbind(as.matrix(sort(rchisq(nsim,1,0))),runif(nsim,0,1)) # chi-square df1 simulations, uniform rvs

Uniftest <- runif(1,0,1)

firstrow <- cbind(LR0,Uniftest)                                       #
use sample LR as LR

NSimM <- rbind(firstrow,NSimM)

Test <- matrix(rep(0,2*(nsim+1)),nrow=(nsim+1))

NSimM <- cbind(NSimM,Test)  

for(i in 2:nsim+1) {                                                 #
indicates the number of simulation above the sample
                                       if (NSimM[i,1]< LR0)NSimM[i,3]<-
1                          # likelihood ratio

         else if(NSimM[i,1]== LR0)if(NSimM[i,2]>= Uniftest)NSimM[i,4]<-1
# if equal, only indicate if rv for simulation
                   }

# is larger that rv for sample LR

Gn <- 1-(sum(NSimM[,3]))/nsim + sum(NSimM[,4])/nsim

pval <- (nsim*Gn+1)/(nsim+1)
#Calculate Monte Carlo p-value

out <- c(pval,confint(fit1))

now <- c(Le,Lw)

LR0

}}  

> test_1 <- dur_ind_test(CDMatrix = Interest,nsim=1000)

Profiling...

> test_1

       A b

42.32406 41.59035

=> likelihood ratio = 1.017641  

Could someone please help?    

http://www.investec.com/EmailDisclaimer/emaildisclaimer.htm

The disclaimer also provides our corporate information and names of our directors as required by law.

The disclaimer is deemed to form part of this message in terms of Section 11 of the Electronic Communications and Transactions Act 25 of 2002. If you cannot access the disclaimer, please obtain a copy thereof from us by sending an email to: disclaimer_at_investec.co.za

        [[alternative HTML version deleted]]



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue 08 Apr 2008 - 10:58:02 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 09 Apr 2008 - 14:30:28 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive