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