dlmodeler.extract: Extract the mean, covariance and prediction intervals for...
In dlmodeler: Generalized Dynamic Linear Modeler

Description Usage Arguments Details Value Author(s) See Also Examples

Extracts the mean (expectation), the variance-covariance matrix, and the prediction intervals for the states and observations of a filtered or smoothed DLM component.

dlmodeler.extract(fs, model, compnames=NULL,
                  type=c("observation","state"),
                  value=c("mean","covariance","interval"),
                  prob=.90)

`fs`	filtered or smoothed `dlmodeler`, as a result from a call to `dlmodeler.filter()` or `dlmodeler.smooth()`.
`model`	object of class `dlmodeler` which was used for filtering or smoothing.
`compnames`	an optional list of components to extract.
`type`	an optional string indicating the type to extract: observation (output, by default) or state.
`value`	an optional string indicating the value to extract: mean (expectation, by default), covariance matrix, or prediction intervals.
`prob`	an optional probability (default = 90%) for the computation of prediction intervals.

A component is a named portion of the state vector matrix which can be extracted with this function. Components are automatically created when DLMs are added together which makes it easier to decompose it later into its building blocks (for example: level+trend+seasonal+cycle).

Let us assume model named m is constructed by adding models named m1 and m2. Typically, m will be constructed with two components named m1 and m1, which can be extracted by this function.

When this function is used with a filtered dlmodeler, it returns the means and covariances of the one-step ahead forecasts for the components:

Zt %*% at = E(y(t) | y(1),y(2)...y(t-1)) for observation means, in the form of a (d,n) matrix.
at = E(alpha(t) | y(1),y(2)...y(t-1)) for state means, in the form of a (m,n) matrix.
Zt %*% Pt %*% t(Zt) + Ht = cov(y(t) | y(1),y(2)...y(t-1)) for observation covariances, in the form of a (d,d,n) array.
Pt = cov(alpha(t) | y(1),y(2)...y(t-1)) for state covariances, in the form of a (m,m,n) array.

When this function is used with a smoothed dlmodeler, it returns the means and covariances of the smoothed components:

Zt %*% at = E(y(t) | y(1),y(2)...y(N)) for observation means, in the form of a (d,n) matrix.
at = E(alpha(t) | y(1),y(2)...y(N)) for state means, in the form of a (m,n) matrix.
Zt %*% Pt %*% t(Zt) + Ht = cov(y(t) | y(1),y(2)...y(N)) for observation covariances, in the form of a (d,d,n) array.
Pt = cov(alpha(t) | y(1),y(2)...y(N)) for state covariances, in the form of a (m,m,n) array.

When the value interval is requested, this function returns a list for each component containing:

mean = the mean (expectaton) for the filtered or smoothed state or observation variable.
lower = lower bound of the prediction interval computed as mean-k*sd, k=-qnorm((1+prob)/2).
upper = upper bound of the prediction interval computed as mean+k*sd, k=-qnorm((1+prob)/2).

Cyrille Szymanski <cnszym@gmail.com>

dlmodeler, dlmodeler.filter, dlmodeler.smooth

## Not run: 
require(dlmodeler)

# generate some data
N <- 365*5
t <- c(1:N,rep(NA,365))
a <- rnorm(N+365,0,.5)
y <- pi + cos(2*pi*t/365.25) + .25*sin(2*pi*t/365.25*3) +
     exp(1)*a + rnorm(N+365,0,.5)

# build a model for this data
m <- dlmodeler.build.polynomial(0,sigmaH=.5,name='level') +
     dlmodeler.build.dseasonal(7,sigmaH=0,name='week') +
     dlmodeler.build.tseasonal(365.25,3,sigmaH=0,name='year') +
     dlmodeler.build.regression(a,sigmaH=0,name='reg')
m$name <- 'mymodel'

system.time(f <- dlmodeler.filter(y, m, raw.result=TRUE))

# extract all the components
m.state.mean <- dlmodeler.extract(f,m,type="state",
                                  value="mean")
m.state.cov <- dlmodeler.extract(f,m,type="state",
                                 value="covariance")
m.obs.mean <- dlmodeler.extract(f,m,type="observation",
                                value="mean")
m.obs.cov <- dlmodeler.extract(f,m,type="observation",
                               value="covariance")
m.obs.int <- dlmodeler.extract(f,m,type="observation",
                               value="interval",prob=.99)

par(mfrow=c(2,1))

# show the one step ahead forecasts & 99% prediction intervals
plot(y,xlim=c(N-10,N+30))
lines(m.obs.int$mymodel$upper[1,],col='light grey')
lines(m.obs.int$mymodel$lower[1,],col='light grey')
lines(m.obs.int$mymodel$mean[1,],col=2)

# see to which values the filter has converged:
m.state.mean$level[,N] # should be close to pi
mean(abs(m.state.mean$week[,N])) # should be close to 0
m.state.mean$year[1,N] # should be close to 1
m.state.mean$year[6,N] # should be close to .25
m.state.mean$reg[,N] # should be close to e

# show the filtered level+year components
plot(m.obs.mean$level[1,]+m.obs.mean$year[1,],
		type='l',ylim=c(pi-2,pi+2),col='light green',
		ylab="smoothed & filtered level+year")

system.time(s <- dlmodeler.smooth(f,m))

# show the smoothed level+year components
s.obs.mean <- dlmodeler.extract(s,m,type="observation",
                                value="mean")
lines(s.obs.mean$level[1,]+s.obs.mean$year[1,],type='l',
		ylim=c(pi-2,pi+2),col='dark green')

## End(Not run)