# Evidently, the model wants W1 to be very small.  It will absorb all variation into
# the weird component.  We could alleviate that problem by putting priors on the
# W's.  But let's try something else instead.  Let's model several series
# simultaneously.  Each series will have its own locally linear component, but they
# will share cyclic components.  That way, the cyclic components should account
# only for features common to the series and the locally linear parts should account
# for what's unique to each series.  Let's look at a few of the series to make sure
# this model looks sensible.

plot.ts ( A535[,2:4] )
plot.ts ( A535[300:400,2:4], plot.type="single" )
# Looks like all three series have the same contribution from heartbeat, but
# they have different signs for respiration

  # First I need a function to build a multidimensional cyclic component.
dlmModTrigMulti <- function ( tau, q=1, a=1, m0, C0, dV=1, dW=0 ) {
  # tau:      the period
  # q:        number of harmonics.  not currently used.
  # n.series: number of time series
  # a:        vector of length n.series of multipliers.  a[j] is the
  #             relative amplitude of series j.
  # m0:       prior mean of theta[0]
  # C0:       prior covariance matrix of theta[0]
  # dV:       vector of the variance of the observational noise
  # dW:       a single number expressing the evolution noise

  mod <- dlmModTrig ( tau=tau, q=q, m0=m0, C0=C0, dV=0, dW=dW )
    # no need for dV here because V will be set later

  mod$FF <- cbind ( a, 0 )
  mod$V  <- diag ( rep(dV, length.out=length(a) ) )

  return ( mod )
} #########################################################
  
  # Try it
dlmModTrigMulti ( 18.75, a=c(3,4,5) )  
dlmModTrigMulti ( 18.75, a=c(3,4,5) ) + dlmModTrigMulti (  2.78, a=c(-2,2,4) ) 
  
  # Also try %+%
dlmModPoly ( order=2, dV=1 ) %+% dlmModPoly ( order=2, dV=2 )
    
  # Build a model for three series.  Each has its own locally linear component.
  # They have the same heartbeat components.  Series #1 has a respiration
  # component opposite to series #2,3.  No weird component.
build.multi1 <- function ( par ) {
	mod <- (  (	    dlmModPoly ( order=2, dV=exp(par[1]), dW=c(0,exp(par[2])) )
                %+% dlmModPoly ( order=2, dV=exp(par[1]), dW=c(0,exp(par[2])) )
                %+% dlmModPoly ( order=2, dV=exp(par[1]), dW=c(0,exp(par[2])) )
              )
              +  dlmModTrigMulti (  2.78, a=c(1,1,1),   dV=0, dW=exp(par[3]) )
              +  dlmModTrigMulti ( 18.75, a=c(1,-1,-1), dV=0, dW=exp(par[4]) )
          )
   return ( mod )
}

  # Try it  
build.multi1 ( 1:4 )


system.time ( mle.multi1 <- dlmMLE ( as.matrix(A535)[,2:4], c(0,0,0,0),
                                     build.multi1 ) )
exp ( mle.multi1$par )
  # fit the model
mod.multi1 <- build.multi1 ( mle.multi1$par )
filt.multi1 <- dlmFilter ( as.matrix(A535)[,2:4], mod.multi1 )
smoo.multi1 <- dlmSmooth ( filt.multi1 )
  # plot the results
par ( mfcol=c(4,3) )
for ( i in 1:3 ) {
	FF <- rep ( 0, 10 )
	FF [ 2*i - 1 ] <- 1
	plot.smoo.component ( smoo.multi1, FF, filt.multi1$y[,i], plot.data=T, xlab="",
	                      ylab="level", main=paste("Series",i), pch=1 )
	FF <- rep ( 0, 10 )
	FF [ 2*i ] <- 1
	plot.smoo.component ( smoo.multi1, FF, plot.data=F, xlab="", ylab="slope" )
	abline ( h=0, lty=3 )
	FF <- rep ( 0, 10 )
	FF [ 7 ] <- 1
	plot.smoo.component ( smoo.multi1, FF, plot.data=F, xlab="", ylab="heart" )
	FF <- rep ( 0, 10 )
	FF [ 9 ] <- 1
	plot.smoo.component ( smoo.multi1, FF, plot.data=F, xlab="", ylab="resp" )
}