1pkgname <- "sspir" 2source(file.path(R.home("share"), "R", "examples-header.R")) 3options(warn = 1) 4options(pager = "console") 5library('sspir') 6 7assign(".oldSearch", search(), pos = 'CheckExEnv') 8cleanEx() 9nameEx("EMalgo") 10### * EMalgo 11 12flush(stderr()); flush(stdout()) 13 14### Name: EMalgo 15### Title: Estimation of variance matrices in a Gaussian state space model 16### Aliases: EMalgo 17### Keywords: models 18 19### ** Examples 20 21## Formulates Gaussian state space model: 22## Trend: local linear 23## Seasonal variation: sine function with frequency = 1 24m1 <- SS( Fmat = function(tt, x, phi) { 25 Fmat <- matrix(NA, nrow=3, ncol=1) 26 Fmat[1,1] <- 1 27 Fmat[2,1] <- cos(2*pi*tt/12) 28 Fmat[3,1] <- sin(2*pi*tt/12) 29 return(Fmat) 30 }, 31 Gmat = function(tt, x, phi) { 32 matrix(c(1,0,0,0,1,0,0,0,1), nrow=3) 33 }, 34 Wmat = matrix(c(0.01,0,0,0,0.1,0,0,0,0.1), nrow=3), 35 Vmat = matrix(1), 36 m0 = matrix(c(0,0,0), nrow=1), 37 C0 = matrix(c(1,0,0,0,0.001,0,0,0,0.001), nrow=3, ncol=3) 38 ) 39 40## Simulates 100 observation from m1 41m1 <- recursion(m1, 100) 42 43## Specifies the correlation structure of W: 44Wstruc <- function(W){ 45 W[1,2:3] <- 0 46 W[2:3,1] <- 0 47 W[2,2] <- (W[2,2]+W[3,3])/2 48 W[3,3] <- W[2,2] 49 W[2,3] <- 0 50 W[3,2] <- 0 51 return(W) 52 } 53 54## Estimstes variances and covariances of W by use of the EM algorithm 55estimates <- EMalgo(m1, Wstruc=Wstruc) 56estimates$ss$Wmat 57 58## Plots estimated model 59plot(estimates$ss) 60 61 62 63cleanEx() 64nameEx("Fkfilter") 65### * Fkfilter 66 67flush(stderr()); flush(stdout()) 68 69### Name: Fkfilter 70### Title: Kalman filtering of Gaussian state space model using'KFS' 71### Aliases: Fkfilter 72### Keywords: models 73 74### ** Examples 75 76## Formulates Gaussian state space model: 77## Trend: local linear 78## Seasonal variation: sine function with frequency = 1 79m1 <- SS( Fmat = function(tt, x, phi) { 80 Fmat <- matrix(NA, nrow=3, ncol=1) 81 Fmat[1,1] <- 1 82 Fmat[2,1] <- cos(2*pi*tt/12) 83 Fmat[3,1] <- sin(2*pi*tt/12) 84 return(Fmat) 85 }, 86 Gmat = function(tt, x, phi) { 87 matrix(c(1,0,0,0,1,0,0,0,1), nrow=3) 88 }, 89 Wmat = function(tt, x, phi) {matrix(c(0.01,0,0,0,0.1,0,0,0,0.1), nrow=3)}, 90 Vmat = function(tt, x, phi) {matrix(1)}, 91 m0 = matrix(c(0,0,0), nrow=1), 92 C0 = matrix(c(1,0,0,0,0.001,0,0,0,0.001), nrow=3, ncol=3) 93 ) 94 95## Simulates 100 observation from m1 96m1 <- recursion(m1, 100) 97 98filtered <- Fkfilter(m1, tvar=c(m1$n, 1, 1, 1)) 99 100plot(filtered$ss) 101 102 103 104 105cleanEx() 106nameEx("Fkfs") 107### * Fkfs 108 109flush(stderr()); flush(stdout()) 110 111### Name: Fkfs 112### Title: Fast Kalman filtering and smoothing of state space models using 113### 'KFS' or 'logLik.SSModel' 114### Aliases: Fkfs 115### Keywords: models 116 117### ** Examples 118 119## Gaussian model to simulate observations from 120m1 <- SS( Fmat = function(tt, x, phi) { 121 Fmat <- matrix(NA, nrow=3, ncol=1) 122 Fmat[1,1] <- 1 123 Fmat[2,1] <- cos(2*pi*tt/12) 124 Fmat[3,1] <- sin(2*pi*tt/12) 125 return(Fmat) 126 }, 127 Gmat = function(tt, x, phi) { 128 matrix(c(1,0,0,0,1,0,0,0,1), nrow=3) 129 }, 130 Wmat = function(tt, x, phi) { 131 matrix(c(0.01,0,0,0,0.1,0,0,0,0.1), nrow=3)}, 132 Vmat = function(tt, x, phi) {matrix(1)}, 133 m0 = matrix(c(0,0,0), nrow=1), 134 C0 = matrix(c(1,0,0,0,0.001,0,0,0,0.001), nrow=3, ncol=3) 135 ) 136 137## Simulates 100 observation from m1 138m1 <- recursion(m1, 100) 139 140## Formulates model og class ssm 141ssm1 <- ssm(m1$y ~ -1 + tvar(polytime(1:m1$n),1) 142 + tvar(polytrig(1:m1$n, 1)), 143 family="gaussian", fit=FALSE) 144 145smoothed <- Fkfs(ssm1$ss, tvar=c(m1$n, 1, 1, 1)) 146 147## Non Gaussian model 148phi.start <- StructTS(log10(UKgas),type="BSM")$coef[c(4,1,2,3)] 149gasmodel <- ssm( log10(UKgas) ~ -1+ tvar(polytime(time,1))+ 150tvar(sumseason(time,4)), phi=phi.start, fit=FALSE) 151 152smoothed <- Fkfs(gasmodel$ss,tvar=c(1,1,1,1)) 153## Trend plot 154## Trend plot 155ts.plot( smoothed$kfas$alphahat[1,] ) 156## Season plot 157ts.plot( smoothed$kfas$alphahat[3,] ) 158 159 160 161cleanEx() 162nameEx("SS") 163### * SS 164 165flush(stderr()); flush(stdout()) 166 167### Name: SS 168### Title: Representation of Gaussian State Space Model 169### Aliases: SS plot.SS print.SS SS-class C0.SS C0<-.SS Fmat.SS Fmat<-.SS 170### Gmat.SS Gmat<-.SS m0.SS m0<-.SS phi.SS phi<-.SS Vmat.SS Vmat<-.SS 171### Wmat.SS Wmat<-.SS 172### Keywords: models 173 174### ** Examples 175 176data(kurit) ## West & Harrison, page 40 177m1 <- SS(y=kurit, 178 Fmat=function(tt,x,phi) return(matrix(1)), 179 Gmat=function(tt,x,phi) return(matrix(1)), 180 Wmat=function(tt,x,phi) return(matrix(5)), ## Alternatively Wmat=matrix(5) 181 Vmat=function(tt,x,phi) return(matrix(100)), ## Alternatively Vmat=matrix(100) 182 m0=matrix(130),C0=matrix(400) 183 ) 184 185plot(m1$y) 186m1.f <- kfilter(m1) 187m1.s <- smoother(m1.f) 188lines(m1.f$m,lty=2,col=2) 189lines(m1.s$m,lty=2,col=2) 190 191## make a model with an intervention at time 10 192m2 <- m1 193Wmat(m2) <- function(tt,x,phi) { 194 if (tt==10) return(matrix(900)) 195 else return(matrix(5)) 196} 197 198m2.f <- kfilter(m2) 199m2.s <- smoother(m2.f) 200lines(m2.f$m,lty=2,col=4) 201lines(m2.s$m,lty=2,col=4) 202 203## Use 'ssm' to construct an SS skeleton 204phi.start <- StructTS(log10(UKgas),type="BSM")$coef[c(4,1,2,3)] 205gasmodel <- ssm( log10(UKgas) ~ -1+ 206 tvar(polytime(time,1))+ 207 tvar(sumseason(time,12)), 208 phi=phi.start) 209 210m0(gasmodel) 211C0(gasmodel) 212phi(gasmodel) 213 214fit <- getFit(gasmodel) 215plot( fit$m[,1:3] ) 216 217 218 219cleanEx() 220nameEx("extended") 221### * extended 222 223flush(stderr()); flush(stdout()) 224 225### Name: extended 226### Title: Iterated Extended Kalman Smoothing 227### Aliases: extended extended.SS print.extended 228### Keywords: models 229 230### ** Examples 231 232data(mumps) 233index <- 1:length(mumps) # use 'index' instead of time 234model <- ssm( mumps ~ -1 + tvar(polytime(index,1)), 235 family=poisson(link=log)) 236results <- getFit(model) 237plot(mumps,type='l',ylab='Number of Cases',xlab='',axes=FALSE) 238lines( exp(results$m[,1]), lwd=2) 239## Alternatives: 240## results2 <- extended(model$ss) 241## results3 <- kfs(model) ## yields the same 242 243 244 245cleanEx() 246nameEx("forecast") 247### * forecast 248 249flush(stderr()); flush(stdout()) 250 251### Name: forecast 252### Title: Forecasts a Gaussian state space model 253### Aliases: forecast 254### Keywords: models 255 256### ** Examples 257 258# Formulates Gaussian state space model 259m1 <- SS( 260 Fmat = function(tt,x,phi){matrix(cos(2*pi*tt/365))}, 261 Gmat = function(tt,x,phi){matrix(1)}, 262 Wmat = function(tt,x,phi){matrix(phi)}, 263 Vmat = function(tt,x,phi){matrix(0.1)}, 264 phi = c(1e-1), 265 m0 = matrix(1), 266 C0 = matrix(100) 267) 268 269# Simulates observations 270set.seed(984375) 271m1 <- recursion(m1, 365) 272plot(m1) 273 274# Fits model 275fit <- kfs(m1) 276plot(fit) 277 278# Change format of variances as these are time-invariant 279fit$Wmat <- matrix(m1$phi) 280fit$Vmat <- matrix(0.1) 281 282# Estimates variances by use of EM algorithm 283est.fit <- EMalgo(fit) 284 285 286# Forecasting 287fcast <- forecast(fit, k=100) 288 289plot(fcast) 290 291 292 293cleanEx() 294nameEx("kfilter") 295### * kfilter 296 297flush(stderr()); flush(stdout()) 298 299### Name: kfilter 300### Title: Kalman filter for Gaussian state space model 301### Aliases: kfilter kfilter.SS filterstep 302### Keywords: models 303 304### ** Examples 305 306data(kurit) 307m1 <- SS(kurit) 308phi(m1) <- c(100,5) 309m0(m1) <- matrix(130) 310C0(m1) <- matrix(400) 311 312m1.f <- kfilter(m1) 313plot(m1$y) 314lines(m1.f$m,lty=2) 315 316 317 318cleanEx() 319nameEx("kfs") 320### * kfs 321 322flush(stderr()); flush(stdout()) 323 324### Name: kfs 325### Title: (Iterated extended) Kalman smoother 326### Aliases: kfs kfs.ssm kfs.SS 327### Keywords: models 328 329### ** Examples 330 331data(mumps) 332index <- 1:length(mumps) 333phi.start <- c(0,0,0.0005,0.0001) 334m3 <- ssm( mumps ~ -1 + tvar(polytime(index,1)) + 335 tvar(polytrig(index,12,1)), 336 family=poisson(link=log), 337 phi=phi.start, C0 = diag(4), 338 fit=FALSE 339) 340 341## The option "fit=FALSE" means that the Kalman Filter/Smoother is not 342## run. 343## At this point you may inspect/change the setup before running 'kfs' 344C0(m3) 345C0(m3) <- 10*diag(4) 346## incorporate possible structural 'jump' at timepoint 10 347Wold <- Wmat(m3) 348Wmat(m3) <- function(tt,x,phi) { 349 W <- Wold(tt,x,phi) 350 if (tt==10) {W[2,2] <- 100*W[2,2]; return(W)} 351 else return(W) 352} 353 354m3.fit <- kfs(m3) 355 356plot(mumps,type='l',ylab='Number of Cases',xlab='') 357lines(exp(m3.fit$m[,1]),type='l',lwd=2) 358 359 360 361cleanEx() 362nameEx("recursion") 363### * recursion 364 365flush(stderr()); flush(stdout()) 366 367### Name: recursion 368### Title: Simulate from a Gaussian state space model 369### Aliases: recursion recursion.SS 370### Keywords: models 371 372### ** Examples 373 374data(kurit) 375m1 <- SS(kurit) 376phi(m1) <- c(100,5) 377m0(m1) <- matrix(130) 378C0(m1) <- matrix(400) 379 380par(mfrow=c(2,1)) 381plot(recursion(m1,100)) 382phi(m1) <- c(5,100) 383plot(recursion(m1,100)) 384 385 386 387graphics::par(get("par.postscript", pos = 'CheckExEnv')) 388cleanEx() 389nameEx("simulate") 390### * simulate 391 392flush(stderr()); flush(stdout()) 393 394### Name: ksimulate 395### Title: Forwards filtering Backwards sampling 396### Aliases: ksimulate ksimulate.SS 397### Keywords: models 398 399### ** Examples 400 401data(kurit) 402m1 <- SS(kurit) 403phi(m1) <- c(100,5) 404m0(m1) <- matrix(130) 405C0(m1) <- matrix(400) 406 407m1 <- kfilter(m1) 408m1.s <- smoother(m1) 409sim <- ksimulate(m1,10) 410 411plot(kurit) 412for (i in 1:10) lines(sim[,,i],lty=2,col=2) 413 414lines(smoother(m1)$m,lwd=2) 415 416 417 418cleanEx() 419nameEx("smoother") 420### * smoother 421 422flush(stderr()); flush(stdout()) 423 424### Name: smoother 425### Title: Kalman smoother for Gaussian state space model 426### Aliases: smoother smoother.SS smootherstep smootherstep.uni 427### print.Smoothed 428### Keywords: models 429 430### ** Examples 431 432data(kurit) 433m1 <- SS(kurit) 434phi(m1) <- c(100,5) 435m0(m1) <- matrix(130) 436C0(m1) <- matrix(400) 437 438m1.s <- smoother(kfilter(m1)) 439plot(m1$y) 440lines(m1.s$m,lty=2) 441 442 443 444cleanEx() 445nameEx("specials") 446### * specials 447 448flush(stderr()); flush(stdout()) 449 450### Name: Specials 451### Title: Special functions used in ssm formulas 452### Aliases: polytrig polytime sumseason season 453### Keywords: models 454 455### ** Examples 456 457polytrig(1:10,degree=2) 458season(1:12,period=4) 459 460 461 462cleanEx() 463nameEx("ssm") 464### * ssm 465 466flush(stderr()); flush(stdout()) 467 468### Name: ssm 469### Title: Define state-space model in a glm-style call. 470### Aliases: ssm kfilter.ssm smoother.ssm C0 C0<- Fmat Fmat<- Gmat Gmat<- 471### m0 m0<- phi phi<- Vmat Vmat<- Wmat Wmat<- C0.ssm C0<-.ssm Fmat.ssm 472### Fmat<-.ssm Gmat.ssm Gmat<-.ssm m0.ssm m0<-.ssm phi.ssm phi<-.ssm 473### Vmat.ssm Vmat<-.ssm Wmat.ssm Wmat<-.ssm getFit print.ssm 474### Keywords: models 475 476### ** Examples 477 478data(vandrivers) 479vandrivers$y <- ts(vandrivers$y,start=1969,frequency=12) 480vd.time <- time(vandrivers$y) 481vd <- ssm( y ~ tvar(1) + seatbelt + sumseason(vd.time,12), 482 family=poisson(link="log"), 483 data=vandrivers, 484 phi = c(1,0.0004), 485 C0=diag(13)*100, 486 fit=FALSE 487 ) 488phi(vd)["(Intercept)"] <- exp(- 2*3.703307 ) 489C0(vd) <- diag(13)*1000 490vd.res <- kfs(vd) 491 492plot( vd.res$m[,1:3] ) 493 494attach(vandrivers) 495plot(y,ylim=c(0,20)) 496lines(exp(vd.res$m[,1]+vd.res$m[,2]*seatbelt),lwd=2 ) 497detach(vandrivers) 498 499 500 501### * <FOOTER> 502### 503cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n") 504grDevices::dev.off() 505### 506### Local variables: *** 507### mode: outline-minor *** 508### outline-regexp: "\\(> \\)?### [*]+" *** 509### End: *** 510quit('no') 511