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