1data(vandrivers) 2#vd.time <- time(vandrivers$y) 3vd <- ssm( y ~ tvar(1) + seatbelt + sumseason(time,12), 4 family=poisson(link="log"), 5 data=vandrivers, 6 phi = c(1,0.0004945505), 7 C0=diag(13)*1000 8 ) 9vd <- getFit(vd) 10 11 12mle <- function(phi.suggest,obj.ext) { 13 cat(":",phi.suggest,":") 14 obj.ext$phi[2] <- phi.suggest 15 ext <- extended(obj.ext) 16 lik <- ext$likelihood 17 cat("loglik(",phi.suggest,") = ",lik) 18 return( lik ) 19} 20 21#res1 <- optim(par=4e-4,fn=mle,obj.ext=vd.res,method=c("L-BFGS-B"), 22# lower=0,upper=Inf, 23# control=list(parscale=1e-4,fnscale=-1), 24# hessian=FALSE) 25 26#res1 <- optim(par=1,fn=mle,obj.ext=vd.res,method=c("SANN"), 27# control=list(fnscale=-1,maxiter=10), 28# hessian=FALSE) 29 30#res1 <- optim(par=vd$ss$phi[2],fn=mle,method=c("Nelder-Mead"), 31# lower=rep(0,1),upper=rep(Inf,1),control=list(fnscale=-1),hessian=FALSE) 32#res1 <- optim(par=vd$ss$phi[2],fn=mle,method=c("Nelder-Mead"), 33# control=list(trace=6,parscale=1e-4,fnscale=-1),hessian=FALSE) 34#vd$ss$phi[2] <- res1$par 35 36 37attach(vandrivers) 38 39#pdf("vandrivers.pdf",width=10,height=6) 40#postscript("vandrivers.pdf",width=10,height=6,horizontal=FALSE) 41par(mfrow=c(1,1)) 42plot(y,ylim=c(0,20),ylab="Vandrivers killed",xlab="Time") 43lines(exp(seatbelt*vd$m[,2] + vd$m[,1]),lwd=2) 44 45vd.sd <- c() 46for (i in 1:length(y)) { 47 thisone <- vd$C[[i]][1:2,1:2] 48 if (seatbelt[i]==0) { vd.sd <- c(vd.sd,sqrt(thisone[1,1])) } 49 else 50 vd.sd <- c(vd.sd,sqrt(sum(thisone))) 51} 52 53lines(exp(seatbelt*vd$m[,2] + vd$m[,1]+2*vd.sd),lty=2) 54lines(exp(seatbelt*vd$m[,2] + vd$m[,1]-2*vd.sd),lty=2) 55#dev.off() 56 57cat("Reduction of casualties:",100*(1-exp(vd$m[1,2])),"%\n") 58detach("vandrivers") 59