dlmodeler.extract(fs, model, compnames=NULL, type=c("observation","state"), value=c("mean","covariance","interval"), prob=.90)dlmodeler, as a result from a call
to dlmodeler.filter() or dlmodeler.smooth().dlmodeler which was used
for filtering or smoothing.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.
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.
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).
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.
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)
Run the code above in your browser using DataLab